Actual source code: matrix.c
1: #define PETSCMAT_DLL
3: /*
4: This is where the abstract matrix operations are defined
5: */
7: #include include/private/matimpl.h
8: #include private/vecimpl.h
10: /* Logging support */
11: PetscCookie MAT_COOKIE = 0;
12: PetscEvent MAT_Mult = 0, MAT_Mults = 0, MAT_MultConstrained = 0, MAT_MultAdd = 0, MAT_MultTranspose = 0;
13: PetscEvent MAT_MultTransposeConstrained = 0, MAT_MultTransposeAdd = 0, MAT_Solve = 0, MAT_Solves = 0, MAT_SolveAdd = 0, MAT_SolveTranspose = 0, MAT_MatSolve = 0;
14: PetscEvent MAT_SolveTransposeAdd = 0, MAT_Relax = 0, MAT_ForwardSolve = 0, MAT_BackwardSolve = 0, MAT_LUFactor = 0, MAT_LUFactorSymbolic = 0;
15: PetscEvent MAT_LUFactorNumeric = 0, MAT_CholeskyFactor = 0, MAT_CholeskyFactorSymbolic = 0, MAT_CholeskyFactorNumeric = 0, MAT_ILUFactor = 0;
16: PetscEvent MAT_ILUFactorSymbolic = 0, MAT_ICCFactorSymbolic = 0, MAT_Copy = 0, MAT_Convert = 0, MAT_Scale = 0, MAT_AssemblyBegin = 0;
17: PetscEvent MAT_AssemblyEnd = 0, MAT_SetValues = 0, MAT_GetValues = 0, MAT_GetRow = 0, MAT_GetSubMatrices = 0, MAT_GetColoring = 0, MAT_GetOrdering = 0;
18: PetscEvent MAT_IncreaseOverlap = 0, MAT_Partitioning = 0, MAT_ZeroEntries = 0, MAT_Load = 0, MAT_View = 0, MAT_AXPY = 0, MAT_FDColoringCreate = 0;
19: PetscEvent MAT_FDColoringApply = 0,MAT_Transpose = 0,MAT_FDColoringFunction = 0;
20: PetscEvent MAT_MatMult = 0, MAT_MatMultSymbolic = 0, MAT_MatMultNumeric = 0;
21: PetscEvent MAT_PtAP = 0, MAT_PtAPSymbolic = 0, MAT_PtAPNumeric = 0;
22: PetscEvent MAT_MatMultTranspose = 0, MAT_MatMultTransposeSymbolic = 0, MAT_MatMultTransposeNumeric = 0;
24: /* nasty global values for MatSetValue() */
25: PetscInt MatSetValue_Row = 0;
26: PetscInt MatSetValue_Column = 0;
27: PetscScalar MatSetValue_Value = 0.0;
31: /*@
32: MatRealPart - Zeros out the imaginary part of the matrix
34: Collective on Mat
36: Input Parameters:
37: . mat - the matrix
39: Level: advanced
42: .seealso: MatImaginaryPart()
43: @*/
45: PetscErrorCode MatRealPart(Mat mat)
46: {
52: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
53: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
54: if (!mat->ops->realpart) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
55: MatPreallocated(mat);
56: (*mat->ops->realpart)(mat);
57: return(0);
58: }
62: /*@
63: MatImaginaryPart - Moves the imaginary part of the matrix to the real part and zeros the imaginary part
65: Collective on Mat
67: Input Parameters:
68: . mat - the matrix
70: Level: advanced
73: .seealso: MatRealPart()
74: @*/
76: PetscErrorCode MatImaginaryPart(Mat mat)
77: {
83: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
84: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
85: if (!mat->ops->imaginarypart) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
86: MatPreallocated(mat);
87: (*mat->ops->imaginarypart)(mat);
88: return(0);
89: }
93: /*@C
94: MatGetRow - Gets a row of a matrix. You MUST call MatRestoreRow()
95: for each row that you get to ensure that your application does
96: not bleed memory.
98: Not Collective
100: Input Parameters:
101: + mat - the matrix
102: - row - the row to get
104: Output Parameters:
105: + ncols - if not NULL, the number of nonzeros in the row
106: . cols - if not NULL, the column numbers
107: - vals - if not NULL, the values
109: Notes:
110: This routine is provided for people who need to have direct access
111: to the structure of a matrix. We hope that we provide enough
112: high-level matrix routines that few users will need it.
114: MatGetRow() always returns 0-based column indices, regardless of
115: whether the internal representation is 0-based (default) or 1-based.
117: For better efficiency, set cols and/or vals to PETSC_NULL if you do
118: not wish to extract these quantities.
120: The user can only examine the values extracted with MatGetRow();
121: the values cannot be altered. To change the matrix entries, one
122: must use MatSetValues().
124: You can only have one call to MatGetRow() outstanding for a particular
125: matrix at a time, per processor. MatGetRow() can only obtain rows
126: associated with the given processor, it cannot get rows from the
127: other processors; for that we suggest using MatGetSubMatrices(), then
128: MatGetRow() on the submatrix. The row indix passed to MatGetRows()
129: is in the global number of rows.
131: Fortran Notes:
132: The calling sequence from Fortran is
133: .vb
134: MatGetRow(matrix,row,ncols,cols,values,ierr)
135: Mat matrix (input)
136: integer row (input)
137: integer ncols (output)
138: integer cols(maxcols) (output)
139: double precision (or double complex) values(maxcols) output
140: .ve
141: where maxcols >= maximum nonzeros in any row of the matrix.
144: Caution:
145: Do not try to change the contents of the output arrays (cols and vals).
146: In some cases, this may corrupt the matrix.
148: Level: advanced
150: Concepts: matrices^row access
152: .seealso: MatRestoreRow(), MatSetValues(), MatGetValues(), MatGetSubmatrices(), MatGetDiagonal()
153: @*/
155: PetscErrorCode MatGetRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
156: {
158: PetscInt incols;
163: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
164: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
165: if (!mat->ops->getrow) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
166: MatPreallocated(mat);
168: (*mat->ops->getrow)(mat,row,&incols,(PetscInt **)cols,(PetscScalar **)vals);
169: if (ncols) *ncols = incols;
171: return(0);
172: }
176: /*@
177: MatConjugate - replaces the matrix values with their complex conjugates
179: Collective on Mat
181: Input Parameters:
182: . mat - the matrix
184: Level: advanced
186: .seealso: VecConjugate()
187: @*/
188: PetscErrorCode MatConjugate(Mat mat)
189: {
194: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
195: if (!mat->ops->conjugate) SETERRQ(PETSC_ERR_SUP,"Not provided for this matrix format, send email to petsc-maint@mcs.anl.gov");
196: (*mat->ops->conjugate)(mat);
197: return(0);
198: }
202: /*@C
203: MatRestoreRow - Frees any temporary space allocated by MatGetRow().
205: Not Collective
207: Input Parameters:
208: + mat - the matrix
209: . row - the row to get
210: . ncols, cols - the number of nonzeros and their columns
211: - vals - if nonzero the column values
213: Notes:
214: This routine should be called after you have finished examining the entries.
216: Fortran Notes:
217: The calling sequence from Fortran is
218: .vb
219: MatRestoreRow(matrix,row,ncols,cols,values,ierr)
220: Mat matrix (input)
221: integer row (input)
222: integer ncols (output)
223: integer cols(maxcols) (output)
224: double precision (or double complex) values(maxcols) output
225: .ve
226: Where maxcols >= maximum nonzeros in any row of the matrix.
228: In Fortran MatRestoreRow() MUST be called after MatGetRow()
229: before another call to MatGetRow() can be made.
231: Level: advanced
233: .seealso: MatGetRow()
234: @*/
235: PetscErrorCode MatRestoreRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
236: {
242: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
243: if (!mat->ops->restorerow) return(0);
244: (*mat->ops->restorerow)(mat,row,ncols,(PetscInt **)cols,(PetscScalar **)vals);
245: return(0);
246: }
250: /*@C
251: MatGetRowUpperTriangular - Sets a flag to enable calls to MatGetRow() for matrix in MATSBAIJ format.
252: You should call MatRestoreRowUpperTriangular() after calling MatGetRow/MatRestoreRow() to disable the flag.
254: Not Collective
256: Input Parameters:
257: + mat - the matrix
259: Notes:
260: The flag is to ensure that users are aware of MatGetRow() only provides the upper trianglular part of the row for the matrices in MATSBAIJ format.
262: Level: advanced
264: Concepts: matrices^row access
266: .seealso: MatRestoreRowRowUpperTriangular()
267: @*/
269: PetscErrorCode MatGetRowUpperTriangular(Mat mat)
270: {
276: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
277: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
278: if (!mat->ops->getrowuppertriangular) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
279: MatPreallocated(mat);
280: (*mat->ops->getrowuppertriangular)(mat);
281: return(0);
282: }
286: /*@C
287: MatRestoreRowUpperTriangular - Disable calls to MatGetRow() for matrix in MATSBAIJ format.
289: Not Collective
291: Input Parameters:
292: + mat - the matrix
294: Notes:
295: This routine should be called after you have finished MatGetRow/MatRestoreRow().
298: Level: advanced
300: .seealso: MatGetRowUpperTriangular()
301: @*/
302: PetscErrorCode MatRestoreRowUpperTriangular(Mat mat)
303: {
308: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
309: if (!mat->ops->restorerowuppertriangular) return(0);
310: (*mat->ops->restorerowuppertriangular)(mat);
311: return(0);
312: }
316: /*@C
317: MatSetOptionsPrefix - Sets the prefix used for searching for all
318: Mat options in the database.
320: Collective on Mat
322: Input Parameter:
323: + A - the Mat context
324: - prefix - the prefix to prepend to all option names
326: Notes:
327: A hyphen (-) must NOT be given at the beginning of the prefix name.
328: The first character of all runtime options is AUTOMATICALLY the hyphen.
330: Level: advanced
332: .keywords: Mat, set, options, prefix, database
334: .seealso: MatSetFromOptions()
335: @*/
336: PetscErrorCode MatSetOptionsPrefix(Mat A,const char prefix[])
337: {
342: PetscObjectSetOptionsPrefix((PetscObject)A,prefix);
343: return(0);
344: }
348: /*@C
349: MatAppendOptionsPrefix - Appends to the prefix used for searching for all
350: Mat options in the database.
352: Collective on Mat
354: Input Parameters:
355: + A - the Mat context
356: - prefix - the prefix to prepend to all option names
358: Notes:
359: A hyphen (-) must NOT be given at the beginning of the prefix name.
360: The first character of all runtime options is AUTOMATICALLY the hyphen.
362: Level: advanced
364: .keywords: Mat, append, options, prefix, database
366: .seealso: MatGetOptionsPrefix()
367: @*/
368: PetscErrorCode MatAppendOptionsPrefix(Mat A,const char prefix[])
369: {
371:
374: PetscObjectAppendOptionsPrefix((PetscObject)A,prefix);
375: return(0);
376: }
380: /*@C
381: MatGetOptionsPrefix - Sets the prefix used for searching for all
382: Mat options in the database.
384: Not Collective
386: Input Parameter:
387: . A - the Mat context
389: Output Parameter:
390: . prefix - pointer to the prefix string used
392: Notes: On the fortran side, the user should pass in a string 'prefix' of
393: sufficient length to hold the prefix.
395: Level: advanced
397: .keywords: Mat, get, options, prefix, database
399: .seealso: MatAppendOptionsPrefix()
400: @*/
401: PetscErrorCode MatGetOptionsPrefix(Mat A,const char *prefix[])
402: {
407: PetscObjectGetOptionsPrefix((PetscObject)A,prefix);
408: return(0);
409: }
413: /*@
414: MatSetUp - Sets up the internal matrix data structures for the later use.
416: Collective on Mat
418: Input Parameters:
419: . A - the Mat context
421: Notes:
422: For basic use of the Mat classes the user need not explicitly call
423: MatSetUp(), since these actions will happen automatically.
425: Level: advanced
427: .keywords: Mat, setup
429: .seealso: MatCreate(), MatDestroy()
430: @*/
431: PetscErrorCode MatSetUp(Mat A)
432: {
437: MatSetUpPreallocation(A);
438: MatSetFromOptions(A);
439: return(0);
440: }
444: /*@C
445: MatView - Visualizes a matrix object.
447: Collective on Mat
449: Input Parameters:
450: + mat - the matrix
451: - viewer - visualization context
453: Notes:
454: The available visualization contexts include
455: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
456: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
457: output where only the first processor opens
458: the file. All other processors send their
459: data to the first processor to print.
460: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
462: The user can open alternative visualization contexts with
463: + PetscViewerASCIIOpen() - Outputs matrix to a specified file
464: . PetscViewerBinaryOpen() - Outputs matrix in binary to a
465: specified file; corresponding input uses MatLoad()
466: . PetscViewerDrawOpen() - Outputs nonzero matrix structure to
467: an X window display
468: - PetscViewerSocketOpen() - Outputs matrix to Socket viewer.
469: Currently only the sequential dense and AIJ
470: matrix types support the Socket viewer.
472: The user can call PetscViewerSetFormat() to specify the output
473: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
474: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
475: + PETSC_VIEWER_ASCII_DEFAULT - default, prints matrix contents
476: . PETSC_VIEWER_ASCII_MATLAB - prints matrix contents in Matlab format
477: . PETSC_VIEWER_ASCII_DENSE - prints entire matrix including zeros
478: . PETSC_VIEWER_ASCII_COMMON - prints matrix contents, using a sparse
479: format common among all matrix types
480: . PETSC_VIEWER_ASCII_IMPL - prints matrix contents, using an implementation-specific
481: format (which is in many cases the same as the default)
482: . PETSC_VIEWER_ASCII_INFO - prints basic information about the matrix
483: size and structure (not the matrix entries)
484: . PETSC_VIEWER_ASCII_INFO_DETAIL - prints more detailed information about
485: the matrix structure
487: Options Database Keys:
488: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
489: . -mat_view_info_detailed - Prints more detailed info
490: . -mat_view - Prints matrix in ASCII format
491: . -mat_view_matlab - Prints matrix in Matlab format
492: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
493: . -display <name> - Sets display name (default is host)
494: . -draw_pause <sec> - Sets number of seconds to pause after display
495: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
496: . -viewer_socket_machine <machine>
497: . -viewer_socket_port <port>
498: . -mat_view_binary - save matrix to file in binary format
499: - -viewer_binary_filename <name>
500: Level: beginner
502: Concepts: matrices^viewing
503: Concepts: matrices^plotting
504: Concepts: matrices^printing
506: .seealso: PetscViewerSetFormat(), PetscViewerASCIIOpen(), PetscViewerDrawOpen(),
507: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), MatLoad()
508: @*/
509: PetscErrorCode MatView(Mat mat,PetscViewer viewer)
510: {
511: PetscErrorCode ierr;
512: PetscInt rows,cols;
513: PetscTruth iascii;
514: const char *cstr;
515: PetscViewerFormat format;
520: if (!viewer) {
521: PetscViewerASCIIGetStdout(mat->comm,&viewer);
522: }
525: if (!mat->assembled) SETERRQ(PETSC_ERR_ORDER,"Must call MatAssemblyBegin/End() before viewing matrix");
526: MatPreallocated(mat);
528: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
529: if (iascii) {
530: PetscViewerGetFormat(viewer,&format);
531: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
532: if (mat->prefix) {
533: PetscViewerASCIIPrintf(viewer,"Matrix Object:(%s)\n",mat->prefix);
534: } else {
535: PetscViewerASCIIPrintf(viewer,"Matrix Object:\n");
536: }
537: PetscViewerASCIIPushTab(viewer);
538: MatGetType(mat,&cstr);
539: MatGetSize(mat,&rows,&cols);
540: PetscViewerASCIIPrintf(viewer,"type=%s, rows=%D, cols=%D\n",cstr,rows,cols);
541: if (mat->ops->getinfo) {
542: MatInfo info;
543: MatGetInfo(mat,MAT_GLOBAL_SUM,&info);
544: PetscViewerASCIIPrintf(viewer,"total: nonzeros=%D, allocated nonzeros=%D\n",
545: (PetscInt)info.nz_used,(PetscInt)info.nz_allocated);
546: }
547: }
548: }
549: if (mat->ops->view) {
550: PetscViewerASCIIPushTab(viewer);
551: (*mat->ops->view)(mat,viewer);
552: PetscViewerASCIIPopTab(viewer);
553: } else if (!iascii) {
554: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
555: }
556: if (iascii) {
557: PetscViewerGetFormat(viewer,&format);
558: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
559: PetscViewerASCIIPopTab(viewer);
560: }
561: }
562: return(0);
563: }
567: /*@C
568: MatScaleSystem - Scale a vector solution and right hand side to
569: match the scaling of a scaled matrix.
570:
571: Collective on Mat
573: Input Parameter:
574: + mat - the matrix
575: . b - right hand side vector (or PETSC_NULL)
576: - x - solution vector (or PETSC_NULL)
579: Notes:
580: For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
581: internally scaled, so this does nothing. For MPIROWBS it
582: permutes and diagonally scales.
584: The KSP methods automatically call this routine when required
585: (via PCPreSolve()) so it is rarely used directly.
587: Level: Developer
589: Concepts: matrices^scaling
591: .seealso: MatUseScaledForm(), MatUnScaleSystem()
592: @*/
593: PetscErrorCode MatScaleSystem(Mat mat,Vec b,Vec x)
594: {
600: MatPreallocated(mat);
604: if (mat->ops->scalesystem) {
605: (*mat->ops->scalesystem)(mat,b,x);
606: }
607: return(0);
608: }
612: /*@C
613: MatUnScaleSystem - Unscales a vector solution and right hand side to
614: match the original scaling of a scaled matrix.
615:
616: Collective on Mat
618: Input Parameter:
619: + mat - the matrix
620: . b - right hand side vector (or PETSC_NULL)
621: - x - solution vector (or PETSC_NULL)
624: Notes:
625: For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
626: internally scaled, so this does nothing. For MPIROWBS it
627: permutes and diagonally scales.
629: The KSP methods automatically call this routine when required
630: (via PCPreSolve()) so it is rarely used directly.
632: Level: Developer
634: .seealso: MatUseScaledForm(), MatScaleSystem()
635: @*/
636: PetscErrorCode MatUnScaleSystem(Mat mat,Vec b,Vec x)
637: {
643: MatPreallocated(mat);
646: if (mat->ops->unscalesystem) {
647: (*mat->ops->unscalesystem)(mat,b,x);
648: }
649: return(0);
650: }
654: /*@C
655: MatUseScaledForm - For matrix storage formats that scale the
656: matrix (for example MPIRowBS matrices are diagonally scaled on
657: assembly) indicates matrix operations (MatMult() etc) are
658: applied using the scaled matrix.
659:
660: Collective on Mat
662: Input Parameter:
663: + mat - the matrix
664: - scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for
665: applying the original matrix
667: Notes:
668: For scaled matrix formats, applying the original, unscaled matrix
669: will be slightly more expensive
671: Level: Developer
673: .seealso: MatScaleSystem(), MatUnScaleSystem()
674: @*/
675: PetscErrorCode MatUseScaledForm(Mat mat,PetscTruth scaled)
676: {
682: MatPreallocated(mat);
683: if (mat->ops->usescaledform) {
684: (*mat->ops->usescaledform)(mat,scaled);
685: }
686: return(0);
687: }
691: /*@
692: MatDestroy - Frees space taken by a matrix.
693:
694: Collective on Mat
696: Input Parameter:
697: . A - the matrix
699: Level: beginner
701: @*/
702: PetscErrorCode MatDestroy(Mat A)
703: {
707: if (--A->refct > 0) return(0);
708: MatPreallocated(A);
709: /* if memory was published with AMS then destroy it */
710: PetscObjectDepublish(A);
711: if (A->ops->destroy) {
712: (*A->ops->destroy)(A);
713: }
714: if (A->mapping) {
715: ISLocalToGlobalMappingDestroy(A->mapping);
716: }
717: if (A->bmapping) {
718: ISLocalToGlobalMappingDestroy(A->bmapping);
719: }
720: PetscFree(A->rmap.range);
721: PetscFree(A->cmap.range);
722: PetscFree(A->spptr);
723: PetscHeaderDestroy(A);
724: return(0);
725: }
729: /*@
730: MatValid - Checks whether a matrix object is valid.
732: Collective on Mat
734: Input Parameter:
735: . m - the matrix to check
737: Output Parameter:
738: flg - flag indicating matrix status, either
739: PETSC_TRUE if matrix is valid, or PETSC_FALSE otherwise.
741: Level: developer
743: Concepts: matrices^validity
744: @*/
745: PetscErrorCode MatValid(Mat m,PetscTruth *flg)
746: {
749: if (!m) *flg = PETSC_FALSE;
750: else if (m->cookie != MAT_COOKIE) *flg = PETSC_FALSE;
751: else *flg = PETSC_TRUE;
752: return(0);
753: }
757: /*@
758: MatSetValues - Inserts or adds a block of values into a matrix.
759: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
760: MUST be called after all calls to MatSetValues() have been completed.
762: Not Collective
764: Input Parameters:
765: + mat - the matrix
766: . v - a logically two-dimensional array of values
767: . m, idxm - the number of rows and their global indices
768: . n, idxn - the number of columns and their global indices
769: - addv - either ADD_VALUES or INSERT_VALUES, where
770: ADD_VALUES adds values to any existing entries, and
771: INSERT_VALUES replaces existing entries with new values
773: Notes:
774: By default the values, v, are row-oriented and unsorted.
775: See MatSetOption() for other options.
777: Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES
778: options cannot be mixed without intervening calls to the assembly
779: routines.
781: MatSetValues() uses 0-based row and column numbers in Fortran
782: as well as in C.
784: Negative indices may be passed in idxm and idxn, these rows and columns are
785: simply ignored. This allows easily inserting element stiffness matrices
786: with homogeneous Dirchlet boundary conditions that you don't want represented
787: in the matrix.
789: Efficiency Alert:
790: The routine MatSetValuesBlocked() may offer much better efficiency
791: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
793: Level: beginner
795: Concepts: matrices^putting entries in
797: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
798: InsertMode, INSERT_VALUES, ADD_VALUES
799: @*/
800: PetscErrorCode MatSetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
801: {
805: if (!m || !n) return(0); /* no values to insert */
811: MatPreallocated(mat);
812: if (mat->insertmode == NOT_SET_VALUES) {
813: mat->insertmode = addv;
814: }
815: #if defined(PETSC_USE_DEBUG)
816: else if (mat->insertmode != addv) {
817: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
818: }
819: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
820: #endif
822: if (mat->assembled) {
823: mat->was_assembled = PETSC_TRUE;
824: mat->assembled = PETSC_FALSE;
825: }
827: if (!mat->ops->setvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
828: (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);
830: return(0);
831: }
836: /*@
837: MatSetValuesRowLocal - Inserts a row (block row for BAIJ matrices) of nonzero
838: values into a matrix
840: Not Collective
842: Input Parameters:
843: + mat - the matrix
844: . row - the (block) row to set
845: - v - a logically two-dimensional array of values
847: Notes:
848: By the values, v, are column-oriented (for the block version) and sorted
850: All the nonzeros in the row must be provided
852: The matrix must have previously had its column indices set
854: The row must belong to this process
856: Level: intermediate
858: Concepts: matrices^putting entries in
860: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
861: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues(), MatSetValuesRow(), MatSetLocalToGlobalMapping()
862: @*/
863: PetscErrorCode MatSetValuesRowLocal(Mat mat,PetscInt row,const PetscScalar v[])
864: {
871: MatSetValuesRow(mat, mat->mapping->indices[row],v);
872: return(0);
873: }
877: /*@
878: MatSetValuesRow - Inserts a row (block row for BAIJ matrices) of nonzero
879: values into a matrix
881: Not Collective
883: Input Parameters:
884: + mat - the matrix
885: . row - the (block) row to set
886: - v - a logically two-dimensional array of values
888: Notes:
889: By the values, v, are column-oriented (for the block version) and sorted
891: All the nonzeros in the row must be provided
893: The matrix must have previously had its column indices set
895: The row must belong to this process
897: Level: intermediate
899: Concepts: matrices^putting entries in
901: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
902: InsertMode, INSERT_VALUES, ADD_VALUES, MatSetValues()
903: @*/
904: PetscErrorCode MatSetValuesRow(Mat mat,PetscInt row,const PetscScalar v[])
905: {
912: #if defined(PETSC_USE_DEBUG)
913: if (mat->insertmode == ADD_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add and insert values");
914: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
915: #endif
916: mat->insertmode = INSERT_VALUES;
918: if (mat->assembled) {
919: mat->was_assembled = PETSC_TRUE;
920: mat->assembled = PETSC_FALSE;
921: }
923: if (!mat->ops->setvaluesrow) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
924: (*mat->ops->setvaluesrow)(mat,row,v);
926: return(0);
927: }
931: /*@
932: MatSetValuesStencil - Inserts or adds a block of values into a matrix.
933: Using structured grid indexing
935: Not Collective
937: Input Parameters:
938: + mat - the matrix
939: . v - a logically two-dimensional array of values
940: . m - number of rows being entered
941: . idxm - grid coordinates (and component number when dof > 1) for matrix rows being entered
942: . n - number of columns being entered
943: . idxn - grid coordinates (and component number when dof > 1) for matrix columns being entered
944: - addv - either ADD_VALUES or INSERT_VALUES, where
945: ADD_VALUES adds values to any existing entries, and
946: INSERT_VALUES replaces existing entries with new values
948: Notes:
949: By default the values, v, are row-oriented and unsorted.
950: See MatSetOption() for other options.
952: Calls to MatSetValuesStencil() with the INSERT_VALUES and ADD_VALUES
953: options cannot be mixed without intervening calls to the assembly
954: routines.
956: The grid coordinates are across the entire grid, not just the local portion
958: MatSetValuesStencil() uses 0-based row and column numbers in Fortran
959: as well as in C.
961: For setting/accessing vector values via array coordinates you can use the DAVecGetArray() routine
963: In order to use this routine you must either obtain the matrix with DAGetMatrix()
964: or call MatSetLocalToGlobalMapping() and MatSetStencil() first.
966: The columns and rows in the stencil passed in MUST be contained within the
967: ghost region of the given process as set with DACreateXXX() or MatSetStencil(). For example,
968: if you create a DA with an overlap of one grid level and on a particular process its first
969: local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
970: first i index you can use in your column and row indices in MatSetStencil() is 5.
972: In Fortran idxm and idxn should be declared as
973: $ MatStencil idxm(4,m),idxn(4,n)
974: and the values inserted using
975: $ idxm(MatStencil_i,1) = i
976: $ idxm(MatStencil_j,1) = j
977: $ idxm(MatStencil_k,1) = k
978: $ idxm(MatStencil_c,1) = c
979: etc
981: For periodic boundary conditions use negative indices for values to the left (below 0; that are to be
982: obtained by wrapping values from right edge). For values to the right of the last entry using that index plus one
983: etc to obtain values that obtained by wrapping the values from the left edge.
985: For indices that don't mean anything for your case (like the k index when working in 2d) or the c index when you have
986: a single value per point) you can skip filling those indices.
988: Inspired by the structured grid interface to the HYPRE package
989: (http://www.llnl.gov/CASC/hypre)
991: Efficiency Alert:
992: The routine MatSetValuesBlockedStencil() may offer much better efficiency
993: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
995: Level: beginner
997: Concepts: matrices^putting entries in
999: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1000: MatSetValues(), MatSetValuesBlockedStencil(), MatSetStencil(), DAGetMatrix(), DAVecGetArray(), MatStencil
1001: @*/
1002: PetscErrorCode MatSetValuesStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1003: {
1005: PetscInt j,i,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1006: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1009: if (!m || !n) return(0); /* no values to insert */
1016: if (m > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 128 rows at a time; trying to set %D",m);
1017: if (n > 256) SETERRQ1(PETSC_ERR_SUP,"Can only set 256 columns at a time; trying to set %D",n);
1019: for (i=0; i<m; i++) {
1020: for (j=0; j<3-sdim; j++) dxm++;
1021: tmp = *dxm++ - starts[0];
1022: for (j=0; j<dim-1; j++) {
1023: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1024: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1025: }
1026: if (mat->stencil.noc) dxm++;
1027: jdxm[i] = tmp;
1028: }
1029: for (i=0; i<n; i++) {
1030: for (j=0; j<3-sdim; j++) dxn++;
1031: tmp = *dxn++ - starts[0];
1032: for (j=0; j<dim-1; j++) {
1033: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1034: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1035: }
1036: if (mat->stencil.noc) dxn++;
1037: jdxn[i] = tmp;
1038: }
1039: MatSetValuesLocal(mat,m,jdxm,n,jdxn,v,addv);
1040: return(0);
1041: }
1045: /*@C
1046: MatSetValuesBlockedStencil - Inserts or adds a block of values into a matrix.
1047: Using structured grid indexing
1049: Not Collective
1051: Input Parameters:
1052: + mat - the matrix
1053: . v - a logically two-dimensional array of values
1054: . m - number of rows being entered
1055: . idxm - grid coordinates for matrix rows being entered
1056: . n - number of columns being entered
1057: . idxn - grid coordinates for matrix columns being entered
1058: - addv - either ADD_VALUES or INSERT_VALUES, where
1059: ADD_VALUES adds values to any existing entries, and
1060: INSERT_VALUES replaces existing entries with new values
1062: Notes:
1063: By default the values, v, are row-oriented and unsorted.
1064: See MatSetOption() for other options.
1066: Calls to MatSetValuesBlockedStencil() with the INSERT_VALUES and ADD_VALUES
1067: options cannot be mixed without intervening calls to the assembly
1068: routines.
1070: The grid coordinates are across the entire grid, not just the local portion
1072: MatSetValuesBlockedStencil() uses 0-based row and column numbers in Fortran
1073: as well as in C.
1075: For setting/accessing vector values via array coordinates you can use the DAVecGetArray() routine
1077: In order to use this routine you must either obtain the matrix with DAGetMatrix()
1078: or call MatSetLocalToGlobalMapping() and MatSetStencil() first.
1080: The columns and rows in the stencil passed in MUST be contained within the
1081: ghost region of the given process as set with DACreateXXX() or MatSetStencil(). For example,
1082: if you create a DA with an overlap of one grid level and on a particular process its first
1083: local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
1084: first i index you can use in your column and row indices in MatSetStencil() is 5.
1086: In Fortran idxm and idxn should be declared as
1087: $ MatStencil idxm(4,m),idxn(4,n)
1088: and the values inserted using
1089: $ idxm(MatStencil_i,1) = i
1090: $ idxm(MatStencil_j,1) = j
1091: $ idxm(MatStencil_k,1) = k
1092: etc
1094: Negative indices may be passed in idxm and idxn, these rows and columns are
1095: simply ignored. This allows easily inserting element stiffness matrices
1096: with homogeneous Dirchlet boundary conditions that you don't want represented
1097: in the matrix.
1099: Inspired by the structured grid interface to the HYPRE package
1100: (http://www.llnl.gov/CASC/hypre)
1102: Level: beginner
1104: Concepts: matrices^putting entries in
1106: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1107: MatSetValues(), MatSetValuesStencil(), MatSetStencil(), DAGetMatrix(), DAVecGetArray(), MatStencil
1108: @*/
1109: PetscErrorCode MatSetValuesBlockedStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
1110: {
1112: PetscInt j,i,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
1113: PetscInt *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);
1116: if (!m || !n) return(0); /* no values to insert */
1123: if (m > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 128 rows at a time; trying to set %D",m);
1124: if (n > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 256 columns at a time; trying to set %D",n);
1126: for (i=0; i<m; i++) {
1127: for (j=0; j<3-sdim; j++) dxm++;
1128: tmp = *dxm++ - starts[0];
1129: for (j=0; j<sdim-1; j++) {
1130: if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1131: else tmp = tmp*dims[j] + *(dxm-1) - starts[j+1];
1132: }
1133: dxm++;
1134: jdxm[i] = tmp;
1135: }
1136: for (i=0; i<n; i++) {
1137: for (j=0; j<3-sdim; j++) dxn++;
1138: tmp = *dxn++ - starts[0];
1139: for (j=0; j<sdim-1; j++) {
1140: if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
1141: else tmp = tmp*dims[j] + *(dxn-1) - starts[j+1];
1142: }
1143: dxn++;
1144: jdxn[i] = tmp;
1145: }
1146: MatSetValuesBlockedLocal(mat,m,jdxm,n,jdxn,v,addv);
1147: return(0);
1148: }
1152: /*@
1153: MatSetStencil - Sets the grid information for setting values into a matrix via
1154: MatSetValuesStencil()
1156: Not Collective
1158: Input Parameters:
1159: + mat - the matrix
1160: . dim - dimension of the grid 1, 2, or 3
1161: . dims - number of grid points in x, y, and z direction, including ghost points on your processor
1162: . starts - starting point of ghost nodes on your processor in x, y, and z direction
1163: - dof - number of degrees of freedom per node
1166: Inspired by the structured grid interface to the HYPRE package
1167: (www.llnl.gov/CASC/hyper)
1169: For matrices generated with DAGetMatrix() this routine is automatically called and so not needed by the
1170: user.
1171:
1172: Level: beginner
1174: Concepts: matrices^putting entries in
1176: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
1177: MatSetValues(), MatSetValuesBlockedStencil(), MatSetValuesStencil()
1178: @*/
1179: PetscErrorCode MatSetStencil(Mat mat,PetscInt dim,const PetscInt dims[],const PetscInt starts[],PetscInt dof)
1180: {
1181: PetscInt i;
1188: mat->stencil.dim = dim + (dof > 1);
1189: for (i=0; i<dim; i++) {
1190: mat->stencil.dims[i] = dims[dim-i-1]; /* copy the values in backwards */
1191: mat->stencil.starts[i] = starts[dim-i-1];
1192: }
1193: mat->stencil.dims[dim] = dof;
1194: mat->stencil.starts[dim] = 0;
1195: mat->stencil.noc = (PetscTruth)(dof == 1);
1196: return(0);
1197: }
1201: /*@
1202: MatSetValuesBlocked - Inserts or adds a block of values into a matrix.
1204: Not Collective
1206: Input Parameters:
1207: + mat - the matrix
1208: . v - a logically two-dimensional array of values
1209: . m, idxm - the number of block rows and their global block indices
1210: . n, idxn - the number of block columns and their global block indices
1211: - addv - either ADD_VALUES or INSERT_VALUES, where
1212: ADD_VALUES adds values to any existing entries, and
1213: INSERT_VALUES replaces existing entries with new values
1215: Notes:
1216: The m and n count the NUMBER of blocks in the row direction and column direction,
1217: NOT the total number of rows/columns; for example, if the block size is 2 and
1218: you are passing in values for rows 2,3,4,5 then m would be 2 (not 4).
1220: By default the values, v, are row-oriented and unsorted. So the layout of
1221: v is the same as for MatSetValues(). See MatSetOption() for other options.
1223: Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
1224: options cannot be mixed without intervening calls to the assembly
1225: routines.
1227: MatSetValuesBlocked() uses 0-based row and column numbers in Fortran
1228: as well as in C.
1230: Negative indices may be passed in idxm and idxn, these rows and columns are
1231: simply ignored. This allows easily inserting element stiffness matrices
1232: with homogeneous Dirchlet boundary conditions that you don't want represented
1233: in the matrix.
1235: Each time an entry is set within a sparse matrix via MatSetValues(),
1236: internal searching must be done to determine where to place the the
1237: data in the matrix storage space. By instead inserting blocks of
1238: entries via MatSetValuesBlocked(), the overhead of matrix assembly is
1239: reduced.
1241: Restrictions:
1242: MatSetValuesBlocked() is currently supported only for the BAIJ and SBAIJ formats
1244: Level: intermediate
1246: Concepts: matrices^putting entries in blocked
1248: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
1249: @*/
1250: PetscErrorCode MatSetValuesBlocked(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
1251: {
1255: if (!m || !n) return(0); /* no values to insert */
1261: MatPreallocated(mat);
1262: if (mat->insertmode == NOT_SET_VALUES) {
1263: mat->insertmode = addv;
1264: }
1265: #if defined(PETSC_USE_DEBUG)
1266: else if (mat->insertmode != addv) {
1267: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1268: }
1269: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1270: #endif
1272: if (mat->assembled) {
1273: mat->was_assembled = PETSC_TRUE;
1274: mat->assembled = PETSC_FALSE;
1275: }
1277: if (!mat->ops->setvaluesblocked) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1278: (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);
1280: return(0);
1281: }
1285: /*@
1286: MatGetValues - Gets a block of values from a matrix.
1288: Not Collective; currently only returns a local block
1290: Input Parameters:
1291: + mat - the matrix
1292: . v - a logically two-dimensional array for storing the values
1293: . m, idxm - the number of rows and their global indices
1294: - n, idxn - the number of columns and their global indices
1296: Notes:
1297: The user must allocate space (m*n PetscScalars) for the values, v.
1298: The values, v, are then returned in a row-oriented format,
1299: analogous to that used by default in MatSetValues().
1301: MatGetValues() uses 0-based row and column numbers in
1302: Fortran as well as in C.
1304: MatGetValues() requires that the matrix has been assembled
1305: with MatAssemblyBegin()/MatAssemblyEnd(). Thus, calls to
1306: MatSetValues() and MatGetValues() CANNOT be made in succession
1307: without intermediate matrix assembly.
1309: Level: advanced
1311: Concepts: matrices^accessing values
1313: .seealso: MatGetRow(), MatGetSubmatrices(), MatSetValues()
1314: @*/
1315: PetscErrorCode MatGetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1316: {
1325: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1326: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1327: if (!mat->ops->getvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1328: MatPreallocated(mat);
1331: (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);
1333: return(0);
1334: }
1338: /*@
1339: MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
1340: the routine MatSetValuesLocal() to allow users to insert matrix entries
1341: using a local (per-processor) numbering.
1343: Not Collective
1345: Input Parameters:
1346: + x - the matrix
1347: - mapping - mapping created with ISLocalToGlobalMappingCreate()
1348: or ISLocalToGlobalMappingCreateIS()
1350: Level: intermediate
1352: Concepts: matrices^local to global mapping
1353: Concepts: local to global mapping^for matrices
1355: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
1356: @*/
1357: PetscErrorCode MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping)
1358: {
1364: if (x->mapping) {
1365: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
1366: }
1367: MatPreallocated(x);
1369: if (x->ops->setlocaltoglobalmapping) {
1370: (*x->ops->setlocaltoglobalmapping)(x,mapping);
1371: } else {
1372: PetscObjectReference((PetscObject)mapping);
1373: x->mapping = mapping;
1374: }
1375: return(0);
1376: }
1380: /*@
1381: MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use
1382: by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
1383: entries using a local (per-processor) numbering.
1385: Not Collective
1387: Input Parameters:
1388: + x - the matrix
1389: - mapping - mapping created with ISLocalToGlobalMappingCreate() or
1390: ISLocalToGlobalMappingCreateIS()
1392: Level: intermediate
1394: Concepts: matrices^local to global mapping blocked
1395: Concepts: local to global mapping^for matrices, blocked
1397: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
1398: MatSetValuesBlocked(), MatSetValuesLocal()
1399: @*/
1400: PetscErrorCode MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping mapping)
1401: {
1407: if (x->bmapping) {
1408: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
1409: }
1410: PetscObjectReference((PetscObject)mapping);
1411: x->bmapping = mapping;
1412: return(0);
1413: }
1417: /*@
1418: MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
1419: using a local ordering of the nodes.
1421: Not Collective
1423: Input Parameters:
1424: + x - the matrix
1425: . nrow, irow - number of rows and their local indices
1426: . ncol, icol - number of columns and their local indices
1427: . y - a logically two-dimensional array of values
1428: - addv - either INSERT_VALUES or ADD_VALUES, where
1429: ADD_VALUES adds values to any existing entries, and
1430: INSERT_VALUES replaces existing entries with new values
1432: Notes:
1433: Before calling MatSetValuesLocal(), the user must first set the
1434: local-to-global mapping by calling MatSetLocalToGlobalMapping().
1436: Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES
1437: options cannot be mixed without intervening calls to the assembly
1438: routines.
1440: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
1441: MUST be called after all calls to MatSetValuesLocal() have been completed.
1443: Level: intermediate
1445: Concepts: matrices^putting entries in with local numbering
1447: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping(),
1448: MatSetValueLocal()
1449: @*/
1450: PetscErrorCode MatSetValuesLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
1451: {
1453: PetscInt irowm[2048],icolm[2048];
1461: MatPreallocated(mat);
1462: if (mat->insertmode == NOT_SET_VALUES) {
1463: mat->insertmode = addv;
1464: }
1465: #if defined(PETSC_USE_DEBUG)
1466: else if (mat->insertmode != addv) {
1467: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1468: }
1469: if (!mat->ops->setvalueslocal && (nrow > 2048 || ncol > 2048)) {
1470: SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %D %D",nrow,ncol);
1471: }
1472: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1473: #endif
1475: if (mat->assembled) {
1476: mat->was_assembled = PETSC_TRUE;
1477: mat->assembled = PETSC_FALSE;
1478: }
1480: if (!mat->ops->setvalueslocal) {
1481: ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm);
1482: ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);
1483: (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);
1484: } else {
1485: (*mat->ops->setvalueslocal)(mat,nrow,irow,ncol,icol,y,addv);
1486: }
1487: mat->same_nonzero = PETSC_FALSE;
1489: return(0);
1490: }
1494: /*@
1495: MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
1496: using a local ordering of the nodes a block at a time.
1498: Not Collective
1500: Input Parameters:
1501: + x - the matrix
1502: . nrow, irow - number of rows and their local indices
1503: . ncol, icol - number of columns and their local indices
1504: . y - a logically two-dimensional array of values
1505: - addv - either INSERT_VALUES or ADD_VALUES, where
1506: ADD_VALUES adds values to any existing entries, and
1507: INSERT_VALUES replaces existing entries with new values
1509: Notes:
1510: Before calling MatSetValuesBlockedLocal(), the user must first set the
1511: local-to-global mapping by calling MatSetLocalToGlobalMappingBlock(),
1512: where the mapping MUST be set for matrix blocks, not for matrix elements.
1514: Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1515: options cannot be mixed without intervening calls to the assembly
1516: routines.
1518: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
1519: MUST be called after all calls to MatSetValuesBlockedLocal() have been completed.
1521: Level: intermediate
1523: Concepts: matrices^putting blocked values in with local numbering
1525: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked()
1526: @*/
1527: PetscErrorCode MatSetValuesBlockedLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
1528: {
1530: PetscInt irowm[2048],icolm[2048];
1538: MatPreallocated(mat);
1539: if (mat->insertmode == NOT_SET_VALUES) {
1540: mat->insertmode = addv;
1541: }
1542: #if defined(PETSC_USE_DEBUG)
1543: else if (mat->insertmode != addv) {
1544: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1545: }
1546: if (!mat->bmapping) {
1547: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with MatSetLocalToGlobalMappingBlock()");
1548: }
1549: if (nrow > 2048 || ncol > 2048) {
1550: SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %D %D",nrow,ncol);
1551: }
1552: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1553: #endif
1555: if (mat->assembled) {
1556: mat->was_assembled = PETSC_TRUE;
1557: mat->assembled = PETSC_FALSE;
1558: }
1560: ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm);
1561: ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm);
1562: (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);
1564: return(0);
1565: }
1567: /* --------------------------------------------------------*/
1570: /*@
1571: MatMult - Computes the matrix-vector product, y = Ax.
1573: Collective on Mat and Vec
1575: Input Parameters:
1576: + mat - the matrix
1577: - x - the vector to be multiplied
1579: Output Parameters:
1580: . y - the result
1582: Notes:
1583: The vectors x and y cannot be the same. I.e., one cannot
1584: call MatMult(A,y,y).
1586: Level: beginner
1588: Concepts: matrix-vector product
1590: .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
1591: @*/
1592: PetscErrorCode MatMult(Mat mat,Vec x,Vec y)
1593: {
1602: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1603: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1604: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1605: #ifndef PETSC_HAVE_CONSTRAINTS
1606: 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);
1607: 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);
1608: 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);
1609: #endif
1610: MatPreallocated(mat);
1612: if (mat->nullsp) {
1613: MatNullSpaceRemove(mat->nullsp,x,&x);
1614: }
1616: if (!mat->ops->mult) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply defined");
1618: (*mat->ops->mult)(mat,x,y);
1621: if (mat->nullsp) {
1622: MatNullSpaceRemove(mat->nullsp,y,PETSC_NULL);
1623: }
1624: PetscObjectStateIncrease((PetscObject)y);
1625: return(0);
1626: }
1630: /*@
1631: MatMultTranspose - Computes matrix transpose times a vector.
1633: Collective on Mat and Vec
1635: Input Parameters:
1636: + mat - the matrix
1637: - x - the vector to be multilplied
1639: Output Parameters:
1640: . y - the result
1642: Notes:
1643: The vectors x and y cannot be the same. I.e., one cannot
1644: call MatMultTranspose(A,y,y).
1646: Level: beginner
1648: Concepts: matrix vector product^transpose
1650: .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd()
1651: @*/
1652: PetscErrorCode MatMultTranspose(Mat mat,Vec x,Vec y)
1653: {
1662: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1663: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1664: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1665: #ifndef PETSC_HAVE_CONSTRAINTS
1666: 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);
1667: 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);
1668: #endif
1669: MatPreallocated(mat);
1671: if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply tranpose defined");
1673: (*mat->ops->multtranspose)(mat,x,y);
1675: PetscObjectStateIncrease((PetscObject)y);
1676: return(0);
1677: }
1681: /*@
1682: MatMultAdd - Computes v3 = v2 + A * v1.
1684: Collective on Mat and Vec
1686: Input Parameters:
1687: + mat - the matrix
1688: - v1, v2 - the vectors
1690: Output Parameters:
1691: . v3 - the result
1693: Notes:
1694: The vectors v1 and v3 cannot be the same. I.e., one cannot
1695: call MatMultAdd(A,v1,v2,v1).
1697: Level: beginner
1699: Concepts: matrix vector product^addition
1701: .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd()
1702: @*/
1703: PetscErrorCode MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1704: {
1714: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1715: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1716: 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);
1717: 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);
1718: 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);
1719: 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);
1720: 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);
1721: if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
1722: MatPreallocated(mat);
1725: (*mat->ops->multadd)(mat,v1,v2,v3);
1727: PetscObjectStateIncrease((PetscObject)v3);
1728: return(0);
1729: }
1733: /*@
1734: MatMultTransposeAdd - Computes v3 = v2 + A' * v1.
1736: Collective on Mat and Vec
1738: Input Parameters:
1739: + mat - the matrix
1740: - v1, v2 - the vectors
1742: Output Parameters:
1743: . v3 - the result
1745: Notes:
1746: The vectors v1 and v3 cannot be the same. I.e., one cannot
1747: call MatMultTransposeAdd(A,v1,v2,v1).
1749: Level: beginner
1751: Concepts: matrix vector product^transpose and addition
1753: .seealso: MatMultTranspose(), MatMultAdd(), MatMult()
1754: @*/
1755: PetscErrorCode MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1756: {
1766: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1767: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1768: if (!mat->ops->multtransposeadd) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1769: if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
1770: 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);
1771: 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);
1772: 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);
1773: MatPreallocated(mat);
1776: (*mat->ops->multtransposeadd)(mat,v1,v2,v3);
1778: PetscObjectStateIncrease((PetscObject)v3);
1779: return(0);
1780: }
1784: /*@
1785: MatMultConstrained - The inner multiplication routine for a
1786: constrained matrix P^T A P.
1788: Collective on Mat and Vec
1790: Input Parameters:
1791: + mat - the matrix
1792: - x - the vector to be multilplied
1794: Output Parameters:
1795: . y - the result
1797: Notes:
1798: The vectors x and y cannot be the same. I.e., one cannot
1799: call MatMult(A,y,y).
1801: Level: beginner
1803: .keywords: matrix, multiply, matrix-vector product, constraint
1804: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
1805: @*/
1806: PetscErrorCode MatMultConstrained(Mat mat,Vec x,Vec y)
1807: {
1814: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1815: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1816: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1817: 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);
1818: 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);
1819: 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);
1822: (*mat->ops->multconstrained)(mat,x,y);
1824: PetscObjectStateIncrease((PetscObject)y);
1826: return(0);
1827: }
1831: /*@
1832: MatMultTransposeConstrained - The inner multiplication routine for a
1833: constrained matrix P^T A^T P.
1835: Collective on Mat and Vec
1837: Input Parameters:
1838: + mat - the matrix
1839: - x - the vector to be multilplied
1841: Output Parameters:
1842: . y - the result
1844: Notes:
1845: The vectors x and y cannot be the same. I.e., one cannot
1846: call MatMult(A,y,y).
1848: Level: beginner
1850: .keywords: matrix, multiply, matrix-vector product, constraint
1851: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
1852: @*/
1853: PetscErrorCode MatMultTransposeConstrained(Mat mat,Vec x,Vec y)
1854: {
1861: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1862: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1863: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1864: 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);
1865: 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);
1868: (*mat->ops->multtransposeconstrained)(mat,x,y);
1870: PetscObjectStateIncrease((PetscObject)y);
1872: return(0);
1873: }
1874: /* ------------------------------------------------------------*/
1877: /*@
1878: MatGetInfo - Returns information about matrix storage (number of
1879: nonzeros, memory, etc.).
1881: Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used
1882: as the flag
1884: Input Parameters:
1885: . mat - the matrix
1887: Output Parameters:
1888: + flag - flag indicating the type of parameters to be returned
1889: (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
1890: MAT_GLOBAL_SUM - sum over all processors)
1891: - info - matrix information context
1893: Notes:
1894: The MatInfo context contains a variety of matrix data, including
1895: number of nonzeros allocated and used, number of mallocs during
1896: matrix assembly, etc. Additional information for factored matrices
1897: is provided (such as the fill ratio, number of mallocs during
1898: factorization, etc.). Much of this info is printed to STDOUT
1899: when using the runtime options
1900: $ -info -mat_view_info
1902: Example for C/C++ Users:
1903: See the file ${PETSC_DIR}/include/petscmat.h for a complete list of
1904: data within the MatInfo context. For example,
1905: .vb
1906: MatInfo info;
1907: Mat A;
1908: double mal, nz_a, nz_u;
1910: MatGetInfo(A,MAT_LOCAL,&info);
1911: mal = info.mallocs;
1912: nz_a = info.nz_allocated;
1913: .ve
1915: Example for Fortran Users:
1916: Fortran users should declare info as a double precision
1917: array of dimension MAT_INFO_SIZE, and then extract the parameters
1918: of interest. See the file ${PETSC_DIR}/include/finclude/petscmat.h
1919: a complete list of parameter names.
1920: .vb
1921: double precision info(MAT_INFO_SIZE)
1922: double precision mal, nz_a
1923: Mat A
1924: integer ierr
1926: call MatGetInfo(A,MAT_LOCAL,info,ierr)
1927: mal = info(MAT_INFO_MALLOCS)
1928: nz_a = info(MAT_INFO_NZ_ALLOCATED)
1929: .ve
1931: Level: intermediate
1933: Concepts: matrices^getting information on
1934:
1935: @*/
1936: PetscErrorCode MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
1937: {
1944: if (!mat->ops->getinfo) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1945: MatPreallocated(mat);
1946: (*mat->ops->getinfo)(mat,flag,info);
1947: return(0);
1948: }
1950: /* ----------------------------------------------------------*/
1953: /*@C
1954: MatILUDTFactor - Performs a drop tolerance ILU factorization.
1956: Collective on Mat
1958: Input Parameters:
1959: + mat - the matrix
1960: . row - row permutation
1961: . col - column permutation
1962: - info - information about the factorization to be done
1964: Output Parameters:
1965: . fact - the factored matrix
1967: Level: developer
1969: Notes:
1970: Most users should employ the simplified KSP interface for linear solvers
1971: instead of working directly with matrix algebra routines such as this.
1972: See, e.g., KSPCreate().
1974: This is currently only supported for the SeqAIJ matrix format using code
1975: from Yousef Saad's SPARSEKIT2 package (translated to C with f2c) and/or
1976: Matlab. SPARSEKIT2 is copyrighted by Yousef Saad with the GNU copyright
1977: and thus can be distributed with PETSc.
1979: Concepts: matrices^ILUDT factorization
1981: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
1982: @*/
1983: PetscErrorCode MatILUDTFactor(Mat mat,IS row,IS col,MatFactorInfo *info,Mat *fact)
1984: {
1994: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1995: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1996: if (!mat->ops->iludtfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1997: MatPreallocated(mat);
1999: (*mat->ops->iludtfactor)(mat,row,col,info,fact);
2001: PetscObjectStateIncrease((PetscObject)*fact);
2003: return(0);
2004: }
2008: /*@
2009: MatLUFactor - Performs in-place LU factorization of matrix.
2011: Collective on Mat
2013: Input Parameters:
2014: + mat - the matrix
2015: . row - row permutation
2016: . col - column permutation
2017: - info - options for factorization, includes
2018: $ fill - expected fill as ratio of original fill.
2019: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2020: $ Run with the option -info to determine an optimal value to use
2022: Notes:
2023: Most users should employ the simplified KSP interface for linear solvers
2024: instead of working directly with matrix algebra routines such as this.
2025: See, e.g., KSPCreate().
2027: This changes the state of the matrix to a factored matrix; it cannot be used
2028: for example with MatSetValues() unless one first calls MatSetUnfactored().
2030: Level: developer
2032: Concepts: matrices^LU factorization
2034: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(),
2035: MatGetOrdering(), MatSetUnfactored(), MatFactorInfo
2037: @*/
2038: PetscErrorCode MatLUFactor(Mat mat,IS row,IS col,MatFactorInfo *info)
2039: {
2048: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2049: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2050: if (!mat->ops->lufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2051: MatPreallocated(mat);
2054: (*mat->ops->lufactor)(mat,row,col,info);
2056: PetscObjectStateIncrease((PetscObject)mat);
2057: return(0);
2058: }
2062: /*@
2063: MatILUFactor - Performs in-place ILU factorization of matrix.
2065: Collective on Mat
2067: Input Parameters:
2068: + mat - the matrix
2069: . row - row permutation
2070: . col - column permutation
2071: - info - structure containing
2072: $ levels - number of levels of fill.
2073: $ expected fill - as ratio of original fill.
2074: $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices
2075: missing diagonal entries)
2077: Notes:
2078: Probably really in-place only when level of fill is zero, otherwise allocates
2079: new space to store factored matrix and deletes previous memory.
2081: Most users should employ the simplified KSP interface for linear solvers
2082: instead of working directly with matrix algebra routines such as this.
2083: See, e.g., KSPCreate().
2085: Level: developer
2087: Concepts: matrices^ILU factorization
2089: .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2090: @*/
2091: PetscErrorCode MatILUFactor(Mat mat,IS row,IS col,MatFactorInfo *info)
2092: {
2101: if (mat->rmap.N != mat->cmap.N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
2102: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2103: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2104: if (!mat->ops->ilufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2105: MatPreallocated(mat);
2108: (*mat->ops->ilufactor)(mat,row,col,info);
2110: PetscObjectStateIncrease((PetscObject)mat);
2111: return(0);
2112: }
2116: /*@
2117: MatLUFactorSymbolic - Performs symbolic LU factorization of matrix.
2118: Call this routine before calling MatLUFactorNumeric().
2120: Collective on Mat
2122: Input Parameters:
2123: + mat - the matrix
2124: . row, col - row and column permutations
2125: - info - options for factorization, includes
2126: $ fill - expected fill as ratio of original fill.
2127: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2128: $ Run with the option -info to determine an optimal value to use
2130: Output Parameter:
2131: . fact - new matrix that has been symbolically factored
2133: Notes:
2134: See the users manual for additional information about
2135: choosing the fill factor for better efficiency.
2137: Most users should employ the simplified KSP interface for linear solvers
2138: instead of working directly with matrix algebra routines such as this.
2139: See, e.g., KSPCreate().
2141: Level: developer
2143: Concepts: matrices^LU symbolic factorization
2145: .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
2146: @*/
2147: PetscErrorCode MatLUFactorSymbolic(Mat mat,IS row,IS col,MatFactorInfo *info,Mat *fact)
2148: {
2158: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2159: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2160: if (!mat->ops->lufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic LU",mat->type_name);
2161: MatPreallocated(mat);
2164: (*mat->ops->lufactorsymbolic)(mat,row,col,info,fact);
2166: PetscObjectStateIncrease((PetscObject)*fact);
2167: return(0);
2168: }
2172: /*@
2173: MatLUFactorNumeric - Performs numeric LU factorization of a matrix.
2174: Call this routine after first calling MatLUFactorSymbolic().
2176: Collective on Mat
2178: Input Parameters:
2179: + mat - the matrix
2180: . info - options for factorization
2181: - fact - the matrix generated for the factor, from MatLUFactorSymbolic()
2183: Notes:
2184: See MatLUFactor() for in-place factorization. See
2185: MatCholeskyFactorNumeric() for the symmetric, positive definite case.
2187: Most users should employ the simplified KSP interface for linear solvers
2188: instead of working directly with matrix algebra routines such as this.
2189: See, e.g., KSPCreate().
2191: Level: developer
2193: Concepts: matrices^LU numeric factorization
2195: .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
2196: @*/
2197: PetscErrorCode MatLUFactorNumeric(Mat mat,MatFactorInfo *info,Mat *fact)
2198: {
2206: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2207: if (mat->rmap.N != (*fact)->rmap.N || mat->cmap.N != (*fact)->cmap.N) {
2208: 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);
2209: }
2210: if (!(*fact)->ops->lufactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2211: MatPreallocated(mat);
2213: (*(*fact)->ops->lufactornumeric)(mat,info,fact);
2216: MatView_Private(*fact);
2217: PetscObjectStateIncrease((PetscObject)*fact);
2218: return(0);
2219: }
2223: /*@
2224: MatCholeskyFactor - Performs in-place Cholesky factorization of a
2225: symmetric matrix.
2227: Collective on Mat
2229: Input Parameters:
2230: + mat - the matrix
2231: . perm - row and column permutations
2232: - f - expected fill as ratio of original fill
2234: Notes:
2235: See MatLUFactor() for the nonsymmetric case. See also
2236: MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().
2238: Most users should employ the simplified KSP interface for linear solvers
2239: instead of working directly with matrix algebra routines such as this.
2240: See, e.g., KSPCreate().
2242: Level: developer
2244: Concepts: matrices^Cholesky factorization
2246: .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
2247: MatGetOrdering()
2249: @*/
2250: PetscErrorCode MatCholeskyFactor(Mat mat,IS perm,MatFactorInfo *info)
2251: {
2259: if (mat->rmap.N != mat->cmap.N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
2260: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2261: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2262: if (!mat->ops->choleskyfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2263: MatPreallocated(mat);
2266: (*mat->ops->choleskyfactor)(mat,perm,info);
2268: PetscObjectStateIncrease((PetscObject)mat);
2269: return(0);
2270: }
2274: /*@
2275: MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
2276: of a symmetric matrix.
2278: Collective on Mat
2280: Input Parameters:
2281: + mat - the matrix
2282: . perm - row and column permutations
2283: - info - options for factorization, includes
2284: $ fill - expected fill as ratio of original fill.
2285: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
2286: $ Run with the option -info to determine an optimal value to use
2288: Output Parameter:
2289: . fact - the factored matrix
2291: Notes:
2292: See MatLUFactorSymbolic() for the nonsymmetric case. See also
2293: MatCholeskyFactor() and MatCholeskyFactorNumeric().
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^Cholesky symbolic factorization
2303: .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
2304: MatGetOrdering()
2306: @*/
2307: PetscErrorCode MatCholeskyFactorSymbolic(Mat mat,IS perm,MatFactorInfo *info,Mat *fact)
2308: {
2317: if (mat->rmap.N != mat->cmap.N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
2318: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2319: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2320: if (!mat->ops->choleskyfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2321: MatPreallocated(mat);
2324: (*mat->ops->choleskyfactorsymbolic)(mat,perm,info,fact);
2326: PetscObjectStateIncrease((PetscObject)*fact);
2327: return(0);
2328: }
2332: /*@
2333: MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
2334: of a symmetric matrix. Call this routine after first calling
2335: MatCholeskyFactorSymbolic().
2337: Collective on Mat
2339: Input Parameter:
2340: . mat - the initial matrix
2341: . info - options for factorization
2342: - fact - the symbolic factor of mat
2344: Output Parameter:
2345: . fact - the factored matrix
2347: Notes:
2348: Most users should employ the simplified KSP interface for linear solvers
2349: instead of working directly with matrix algebra routines such as this.
2350: See, e.g., KSPCreate().
2352: Level: developer
2354: Concepts: matrices^Cholesky numeric factorization
2356: .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
2357: @*/
2358: PetscErrorCode MatCholeskyFactorNumeric(Mat mat,MatFactorInfo *info,Mat *fact)
2359: {
2367: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2368: if (!(*fact)->ops->choleskyfactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2369: if (mat->rmap.N != (*fact)->rmap.N || mat->cmap.N != (*fact)->cmap.N) {
2370: 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);
2371: }
2372: MatPreallocated(mat);
2375: (*(*fact)->ops->choleskyfactornumeric)(mat,info,fact);
2378: MatView_Private(*fact);
2379: PetscObjectStateIncrease((PetscObject)*fact);
2380: return(0);
2381: }
2383: /* ----------------------------------------------------------------*/
2386: /*@
2387: MatSolve - Solves A x = b, given a factored matrix.
2389: Collective on Mat and Vec
2391: Input Parameters:
2392: + mat - the factored matrix
2393: - b - the right-hand-side vector
2395: Output Parameter:
2396: . x - the result vector
2398: Notes:
2399: The vectors b and x cannot be the same. I.e., one cannot
2400: call MatSolve(A,x,x).
2402: Notes:
2403: Most users should employ the simplified KSP interface for linear solvers
2404: instead of working directly with matrix algebra routines such as this.
2405: See, e.g., KSPCreate().
2407: Level: developer
2409: Concepts: matrices^triangular solves
2411: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd()
2412: @*/
2413: PetscErrorCode MatSolve(Mat mat,Vec b,Vec x)
2414: {
2424: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2425: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2426: 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);
2427: 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);
2428: 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);
2429: if (!mat->rmap.N && !mat->cmap.N) return(0);
2430: if (!mat->ops->solve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2431: MatPreallocated(mat);
2434: (*mat->ops->solve)(mat,b,x);
2436: PetscObjectStateIncrease((PetscObject)x);
2437: return(0);
2438: }
2442: /*@
2443: MatMatSolve - Solves A X = B, given a factored matrix.
2445: Collective on Mat
2447: Input Parameters:
2448: + mat - the factored matrix
2449: - b - the right-hand-side vector
2451: Output Parameter:
2452: . x - the result vector
2454: Notes:
2455: The vectors b and x cannot be the same. I.e., one cannot
2456: call MatMatSolve(A,x,x).
2458: Notes:
2459: Most users should employ the simplified KSP interface for linear solvers
2460: instead of working directly with matrix algebra routines such as this.
2461: See, e.g., KSPCreate().
2463: Level: developer
2465: Concepts: matrices^triangular solves
2467: .seealso: MatMatSolveAdd(), MatMatSolveTranspose(), MatMatSolveTransposeAdd()
2468: @*/
2469: PetscErrorCode MatMatSolve(Mat A,Mat B,Mat X)
2470: {
2480: if (X == B) SETERRQ(PETSC_ERR_ARG_IDN,"X and B must be different matrices");
2481: if (!A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2482: 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);
2483: 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);
2484: 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);
2485: if (!A->rmap.N && !A->cmap.N) return(0);
2486: if (!A->ops->matsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",A->type_name);
2487: MatPreallocated(A);
2490: (*A->ops->matsolve)(A,B,X);
2492: PetscObjectStateIncrease((PetscObject)X);
2493: return(0);
2494: }
2499: /* @
2500: MatForwardSolve - Solves L x = b, given a factored matrix, A = LU, or
2501: U^T*D^(1/2) x = b, given a factored symmetric matrix, A = U^T*D*U,
2503: Collective on Mat and Vec
2505: Input Parameters:
2506: + mat - the factored matrix
2507: - b - the right-hand-side vector
2509: Output Parameter:
2510: . x - the result vector
2512: Notes:
2513: MatSolve() should be used for most applications, as it performs
2514: a forward solve followed by a backward solve.
2516: The vectors b and x cannot be the same, i.e., one cannot
2517: call MatForwardSolve(A,x,x).
2519: For matrix in seqsbaij format with block size larger than 1,
2520: the diagonal blocks are not implemented as D = D^(1/2) * D^(1/2) yet.
2521: MatForwardSolve() solves U^T*D y = b, and
2522: MatBackwardSolve() solves U x = y.
2523: Thus they do not provide a symmetric preconditioner.
2525: Most users should employ the simplified KSP interface for linear solvers
2526: instead of working directly with matrix algebra routines such as this.
2527: See, e.g., KSPCreate().
2529: Level: developer
2531: Concepts: matrices^forward solves
2533: .seealso: MatSolve(), MatBackwardSolve()
2534: @ */
2535: PetscErrorCode MatForwardSolve(Mat mat,Vec b,Vec x)
2536: {
2546: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2547: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2548: if (!mat->ops->forwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2549: 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);
2550: 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);
2551: 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);
2552: MatPreallocated(mat);
2554: (*mat->ops->forwardsolve)(mat,b,x);
2556: PetscObjectStateIncrease((PetscObject)x);
2557: return(0);
2558: }
2562: /* @
2563: MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.
2564: D^(1/2) U x = b, given a factored symmetric matrix, A = U^T*D*U,
2566: Collective on Mat and Vec
2568: Input Parameters:
2569: + mat - the factored matrix
2570: - b - the right-hand-side vector
2572: Output Parameter:
2573: . x - the result vector
2575: Notes:
2576: MatSolve() should be used for most applications, as it performs
2577: a forward solve followed by a backward solve.
2579: The vectors b and x cannot be the same. I.e., one cannot
2580: call MatBackwardSolve(A,x,x).
2582: For matrix in seqsbaij format with block size larger than 1,
2583: the diagonal blocks are not implemented as D = D^(1/2) * D^(1/2) yet.
2584: MatForwardSolve() solves U^T*D y = b, and
2585: MatBackwardSolve() solves U x = y.
2586: Thus they do not provide a symmetric preconditioner.
2588: Most users should employ the simplified KSP interface for linear solvers
2589: instead of working directly with matrix algebra routines such as this.
2590: See, e.g., KSPCreate().
2592: Level: developer
2594: Concepts: matrices^backward solves
2596: .seealso: MatSolve(), MatForwardSolve()
2597: @ */
2598: PetscErrorCode MatBackwardSolve(Mat mat,Vec b,Vec x)
2599: {
2609: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2610: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2611: if (!mat->ops->backwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2612: 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);
2613: 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);
2614: 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);
2615: MatPreallocated(mat);
2618: (*mat->ops->backwardsolve)(mat,b,x);
2620: PetscObjectStateIncrease((PetscObject)x);
2621: return(0);
2622: }
2626: /*@
2627: MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.
2629: Collective on Mat and Vec
2631: Input Parameters:
2632: + mat - the factored matrix
2633: . b - the right-hand-side vector
2634: - y - the vector to be added to
2636: Output Parameter:
2637: . x - the result vector
2639: Notes:
2640: The vectors b and x cannot be the same. I.e., one cannot
2641: call MatSolveAdd(A,x,y,x).
2643: Most users should employ the simplified KSP interface for linear solvers
2644: instead of working directly with matrix algebra routines such as this.
2645: See, e.g., KSPCreate().
2647: Level: developer
2649: Concepts: matrices^triangular solves
2651: .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd()
2652: @*/
2653: PetscErrorCode MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
2654: {
2655: PetscScalar one = 1.0;
2656: Vec tmp;
2668: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2669: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2670: 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);
2671: 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);
2672: 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);
2673: 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);
2674: 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);
2675: MatPreallocated(mat);
2678: if (mat->ops->solveadd) {
2679: (*mat->ops->solveadd)(mat,b,y,x);
2680: } else {
2681: /* do the solve then the add manually */
2682: if (x != y) {
2683: MatSolve(mat,b,x);
2684: VecAXPY(x,one,y);
2685: } else {
2686: VecDuplicate(x,&tmp);
2687: PetscLogObjectParent(mat,tmp);
2688: VecCopy(x,tmp);
2689: MatSolve(mat,b,x);
2690: VecAXPY(x,one,tmp);
2691: VecDestroy(tmp);
2692: }
2693: }
2695: PetscObjectStateIncrease((PetscObject)x);
2696: return(0);
2697: }
2701: /*@
2702: MatSolveTranspose - Solves A' x = b, given a factored matrix.
2704: Collective on Mat and Vec
2706: Input Parameters:
2707: + mat - the factored matrix
2708: - b - the right-hand-side vector
2710: Output Parameter:
2711: . x - the result vector
2713: Notes:
2714: The vectors b and x cannot be the same. I.e., one cannot
2715: call MatSolveTranspose(A,x,x).
2717: Most users should employ the simplified KSP interface for linear solvers
2718: instead of working directly with matrix algebra routines such as this.
2719: See, e.g., KSPCreate().
2721: Level: developer
2723: Concepts: matrices^triangular solves
2725: .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd()
2726: @*/
2727: PetscErrorCode MatSolveTranspose(Mat mat,Vec b,Vec x)
2728: {
2738: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2739: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2740: if (!mat->ops->solvetranspose) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s",mat->type_name);
2741: 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);
2742: 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);
2743: MatPreallocated(mat);
2745: (*mat->ops->solvetranspose)(mat,b,x);
2747: PetscObjectStateIncrease((PetscObject)x);
2748: return(0);
2749: }
2753: /*@
2754: MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a
2755: factored matrix.
2757: Collective on Mat and Vec
2759: Input Parameters:
2760: + mat - the factored matrix
2761: . b - the right-hand-side vector
2762: - y - the vector to be added to
2764: Output Parameter:
2765: . x - the result vector
2767: Notes:
2768: The vectors b and x cannot be the same. I.e., one cannot
2769: call MatSolveTransposeAdd(A,x,y,x).
2771: Most users should employ the simplified KSP interface for linear solvers
2772: instead of working directly with matrix algebra routines such as this.
2773: See, e.g., KSPCreate().
2775: Level: developer
2777: Concepts: matrices^triangular solves
2779: .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose()
2780: @*/
2781: PetscErrorCode MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x)
2782: {
2783: PetscScalar one = 1.0;
2785: Vec tmp;
2796: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2797: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2798: 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);
2799: 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);
2800: 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);
2801: 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);
2802: MatPreallocated(mat);
2805: if (mat->ops->solvetransposeadd) {
2806: (*mat->ops->solvetransposeadd)(mat,b,y,x);
2807: } else {
2808: /* do the solve then the add manually */
2809: if (x != y) {
2810: MatSolveTranspose(mat,b,x);
2811: VecAXPY(x,one,y);
2812: } else {
2813: VecDuplicate(x,&tmp);
2814: PetscLogObjectParent(mat,tmp);
2815: VecCopy(x,tmp);
2816: MatSolveTranspose(mat,b,x);
2817: VecAXPY(x,one,tmp);
2818: VecDestroy(tmp);
2819: }
2820: }
2822: PetscObjectStateIncrease((PetscObject)x);
2823: return(0);
2824: }
2825: /* ----------------------------------------------------------------*/
2829: /*@
2830: MatRelax - Computes relaxation (SOR, Gauss-Seidel) sweeps.
2832: Collective on Mat and Vec
2834: Input Parameters:
2835: + mat - the matrix
2836: . b - the right hand side
2837: . omega - the relaxation factor
2838: . flag - flag indicating the type of SOR (see below)
2839: . shift - diagonal shift
2840: - its - the number of iterations
2841: - lits - the number of local iterations
2843: Output Parameters:
2844: . x - the solution (can contain an initial guess)
2846: SOR Flags:
2847: . SOR_FORWARD_SWEEP - forward SOR
2848: . SOR_BACKWARD_SWEEP - backward SOR
2849: . SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
2850: . SOR_LOCAL_FORWARD_SWEEP - local forward SOR
2851: . SOR_LOCAL_BACKWARD_SWEEP - local forward SOR
2852: . SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
2853: . SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies
2854: upper/lower triangular part of matrix to
2855: vector (with omega)
2856: . SOR_ZERO_INITIAL_GUESS - zero initial guess
2858: Notes:
2859: SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
2860: SOR_LOCAL_SYMMETRIC_SWEEP perform separate independent smoothings
2861: on each processor.
2863: Application programmers will not generally use MatRelax() directly,
2864: but instead will employ the KSP/PC interface.
2866: Notes for Advanced Users:
2867: The flags are implemented as bitwise inclusive or operations.
2868: For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
2869: to specify a zero initial guess for SSOR.
2871: Most users should employ the simplified KSP interface for linear solvers
2872: instead of working directly with matrix algebra routines such as this.
2873: See, e.g., KSPCreate().
2875: See also, MatPBRelax(). This routine will automatically call the point block
2876: version if the point version is not available.
2878: Level: developer
2880: Concepts: matrices^relaxation
2881: Concepts: matrices^SOR
2882: Concepts: matrices^Gauss-Seidel
2884: @*/
2885: PetscErrorCode MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec x)
2886: {
2896: if (!mat->ops->relax && !mat->ops->pbrelax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2897: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2898: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2899: 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);
2900: 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);
2901: 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);
2902: MatPreallocated(mat);
2904: if (mat->ops->relax) {
2905: ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,lits,x);
2906: } else {
2907: ierr =(*mat->ops->pbrelax)(mat,b,omega,flag,shift,its,lits,x);
2908: }
2910: PetscObjectStateIncrease((PetscObject)x);
2911: return(0);
2912: }
2916: /*@
2917: MatPBRelax - Computes relaxation (SOR, Gauss-Seidel) sweeps.
2919: Collective on Mat and Vec
2921: See MatRelax() for usage
2923: For multi-component PDEs where the Jacobian is stored in a point block format
2924: (with the PETSc BAIJ matrix formats) the relaxation is done one point block at
2925: a time. That is, the small (for example, 4 by 4) blocks along the diagonal are solved
2926: simultaneously (that is a 4 by 4 linear solve is done) to update all the values at a point.
2928: Level: developer
2930: @*/
2931: PetscErrorCode MatPBRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec x)
2932: {
2942: if (!mat->ops->pbrelax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2943: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2944: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2945: 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);
2946: 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);
2947: 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);
2948: MatPreallocated(mat);
2951: ierr =(*mat->ops->pbrelax)(mat,b,omega,flag,shift,its,lits,x);
2953: PetscObjectStateIncrease((PetscObject)x);
2954: return(0);
2955: }
2959: /*
2960: Default matrix copy routine.
2961: */
2962: PetscErrorCode MatCopy_Basic(Mat A,Mat B,MatStructure str)
2963: {
2964: PetscErrorCode ierr;
2965: PetscInt i,rstart,rend,nz;
2966: const PetscInt *cwork;
2967: const PetscScalar *vwork;
2970: if (B->assembled) {
2971: MatZeroEntries(B);
2972: }
2973: MatGetOwnershipRange(A,&rstart,&rend);
2974: for (i=rstart; i<rend; i++) {
2975: MatGetRow(A,i,&nz,&cwork,&vwork);
2976: MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);
2977: MatRestoreRow(A,i,&nz,&cwork,&vwork);
2978: }
2979: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2980: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2981: PetscObjectStateIncrease((PetscObject)B);
2982: return(0);
2983: }
2987: /*@
2988: MatCopy - Copys a matrix to another matrix.
2990: Collective on Mat
2992: Input Parameters:
2993: + A - the matrix
2994: - str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN
2996: Output Parameter:
2997: . B - where the copy is put
2999: Notes:
3000: If you use SAME_NONZERO_PATTERN then the two matrices had better have the
3001: same nonzero pattern or the routine will crash.
3003: MatCopy() copies the matrix entries of a matrix to another existing
3004: matrix (after first zeroing the second matrix). A related routine is
3005: MatConvert(), which first creates a new matrix and then copies the data.
3007: Level: intermediate
3008:
3009: Concepts: matrices^copying
3011: .seealso: MatConvert(), MatDuplicate()
3013: @*/
3014: PetscErrorCode MatCopy(Mat A,Mat B,MatStructure str)
3015: {
3023: MatPreallocated(B);
3025: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3026: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3027: 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);
3028: MatPreallocated(A);
3031: if (A->ops->copy) {
3032: (*A->ops->copy)(A,B,str);
3033: } else { /* generic conversion */
3034: MatCopy_Basic(A,B,str);
3035: }
3036: if (A->mapping) {
3037: if (B->mapping) {ISLocalToGlobalMappingDestroy(B->mapping);B->mapping = 0;}
3038: MatSetLocalToGlobalMapping(B,A->mapping);
3039: }
3040: if (A->bmapping) {
3041: if (B->bmapping) {ISLocalToGlobalMappingDestroy(B->bmapping);B->bmapping = 0;}
3042: MatSetLocalToGlobalMappingBlock(B,A->mapping);
3043: }
3045: PetscObjectStateIncrease((PetscObject)B);
3046: return(0);
3047: }
3051: /*@C
3052: MatConvert - Converts a matrix to another matrix, either of the same
3053: or different type.
3055: Collective on Mat
3057: Input Parameters:
3058: + mat - the matrix
3059: . newtype - new matrix type. Use MATSAME to create a new matrix of the
3060: same type as the original matrix.
3061: - reuse - denotes if the destination matrix is to be created or reused. Currently
3062: MAT_REUSE_MATRIX is only supported for inplace conversion, otherwise use
3063: MAT_INITIAL_MATRIX.
3064: Output Parameter:
3065: . M - pointer to place new matrix
3067: Notes:
3068: MatConvert() first creates a new matrix and then copies the data from
3069: the first matrix. A related routine is MatCopy(), which copies the matrix
3070: entries of one matrix to another already existing matrix context.
3072: Level: intermediate
3074: Concepts: matrices^converting between storage formats
3076: .seealso: MatCopy(), MatDuplicate()
3077: @*/
3078: PetscErrorCode MatConvert(Mat mat, MatType newtype,MatReuse reuse,Mat *M)
3079: {
3080: PetscErrorCode ierr;
3081: PetscTruth sametype,issame,flg;
3082: char convname[256],mtype[256];
3083: Mat B;
3089: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3090: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3091: MatPreallocated(mat);
3093: PetscOptionsGetString(PETSC_NULL,"-matconvert_type",mtype,256,&flg);
3094: if (flg) {
3095: newtype = mtype;
3096: }
3098:
3099: PetscTypeCompare((PetscObject)mat,newtype,&sametype);
3100: PetscStrcmp(newtype,"same",&issame);
3101: if ((reuse==MAT_REUSE_MATRIX) && (mat != *M)) {
3102: SETERRQ(PETSC_ERR_SUP,"MAT_REUSE_MATRIX only supported for in-place conversion currently");
3103: }
3104: if ((sametype || issame) && (reuse==MAT_INITIAL_MATRIX) && mat->ops->duplicate) {
3105: (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);
3106: } else {
3107: PetscErrorCode (*conv)(Mat, MatType,MatReuse,Mat*)=PETSC_NULL;
3108: const char *prefix[3] = {"seq","mpi",""};
3109: PetscInt i;
3111: /*
3112: Order of precedence:
3113: 1) See if a specialized converter is known to the current matrix.
3114: 2) See if a specialized converter is known to the desired matrix class.
3115: 3) See if a good general converter is registered for the desired class
3116: (as of 6/27/03 only MATMPIADJ falls into this category).
3117: 4) See if a good general converter is known for the current matrix.
3118: 5) Use a really basic converter.
3119: */
3120: for (i=0; i<3; i++) {
3121: /* 1) See if a specialized converter is known to the current matrix and the desired class */
3122: PetscStrcpy(convname,"MatConvert_");
3123: PetscStrcat(convname,mat->type_name);
3124: PetscStrcat(convname,"_");
3125: PetscStrcat(convname,prefix[i]);
3126: PetscStrcat(convname,newtype);
3127: PetscStrcat(convname,"_C");
3128: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
3129: if (conv) goto foundconv;
3130: }
3132: /* 2) See if a specialized converter is known to the desired matrix class. */
3133: MatCreate(mat->comm,&B);
3134: MatSetSizes(B,mat->rmap.n,mat->cmap.n,mat->rmap.N,mat->cmap.N);
3135: MatSetType(B,newtype);
3136: for (i=0; i<3; i++) {
3137: PetscStrcpy(convname,"MatConvert_");
3138: PetscStrcat(convname,mat->type_name);
3139: PetscStrcat(convname,"_");
3140: PetscStrcat(convname,prefix[i]);
3141: PetscStrcat(convname,newtype);
3142: PetscStrcat(convname,"_C");
3143: PetscObjectQueryFunction((PetscObject)B,convname,(void (**)(void))&conv);
3144: if (conv) {
3145: MatDestroy(B);
3146: goto foundconv;
3147: }
3148: }
3150: /* 3) See if a good general converter is registered for the desired class */
3151: conv = B->ops->convertfrom;
3152: MatDestroy(B);
3153: if (conv) goto foundconv;
3155: /* 4) See if a good general converter is known for the current matrix */
3156: if (mat->ops->convert) {
3157: conv = mat->ops->convert;
3158: }
3159: if (conv) goto foundconv;
3161: /* 5) Use a really basic converter. */
3162: conv = MatConvert_Basic;
3164: foundconv:
3165: (*conv)(mat,newtype,reuse,M);
3166: }
3167: B = *M;
3169: PetscObjectStateIncrease((PetscObject)B);
3170: return(0);
3171: }
3176: /*@
3177: MatDuplicate - Duplicates a matrix including the non-zero structure.
3179: Collective on Mat
3181: Input Parameters:
3182: + mat - the matrix
3183: - op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero
3184: values as well or not
3186: Output Parameter:
3187: . M - pointer to place new matrix
3189: Level: intermediate
3191: Concepts: matrices^duplicating
3193: .seealso: MatCopy(), MatConvert()
3194: @*/
3195: PetscErrorCode MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
3196: {
3198: Mat B;
3204: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3205: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3206: MatPreallocated(mat);
3208: *M = 0;
3210: if (!mat->ops->duplicate) {
3211: SETERRQ(PETSC_ERR_SUP,"Not written for this matrix type");
3212: }
3213: (*mat->ops->duplicate)(mat,op,M);
3214: B = *M;
3215: if (mat->mapping) {
3216: MatSetLocalToGlobalMapping(B,mat->mapping);
3217: }
3218: if (mat->bmapping) {
3219: MatSetLocalToGlobalMappingBlock(B,mat->bmapping);
3220: }
3221: PetscMapCopy(mat->comm,&mat->rmap,&B->rmap);
3222: PetscMapCopy(mat->comm,&mat->cmap,&B->cmap);
3223:
3225: PetscObjectStateIncrease((PetscObject)B);
3226: return(0);
3227: }
3231: /*@
3232: MatGetDiagonal - Gets the diagonal of a matrix.
3234: Collective on Mat and Vec
3236: Input Parameters:
3237: + mat - the matrix
3238: - v - the vector for storing the diagonal
3240: Output Parameter:
3241: . v - the diagonal of the matrix
3243: Notes:
3244: For the SeqAIJ matrix format, this routine may also be called
3245: on a LU factored matrix; in that case it routines the reciprocal of
3246: the diagonal entries in U. It returns the entries permuted by the
3247: row and column permutation used during the symbolic factorization.
3249: Level: intermediate
3251: Concepts: matrices^accessing diagonals
3253: .seealso: MatGetRow(), MatGetSubmatrices(), MatGetSubmatrix(), MatGetRowMax()
3254: @*/
3255: PetscErrorCode MatGetDiagonal(Mat mat,Vec v)
3256: {
3263: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3264: if (!mat->ops->getdiagonal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3265: MatPreallocated(mat);
3267: (*mat->ops->getdiagonal)(mat,v);
3268: PetscObjectStateIncrease((PetscObject)v);
3269: return(0);
3270: }
3274: /*@
3275: MatGetRowMax - Gets the maximum value (in absolute value) of each
3276: row of the matrix
3278: Collective on Mat and Vec
3280: Input Parameters:
3281: . mat - the matrix
3283: Output Parameter:
3284: . v - the vector for storing the maximums
3286: Level: intermediate
3288: Concepts: matrices^getting row maximums
3290: .seealso: MatGetDiagonal(), MatGetSubmatrices(), MatGetSubmatrix()
3291: @*/
3292: PetscErrorCode MatGetRowMax(Mat mat,Vec v)
3293: {
3300: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3301: if (!mat->ops->getrowmax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3302: MatPreallocated(mat);
3304: (*mat->ops->getrowmax)(mat,v);
3305: PetscObjectStateIncrease((PetscObject)v);
3306: return(0);
3307: }
3311: /*@C
3312: MatTranspose - Computes an in-place or out-of-place transpose of a matrix.
3314: Collective on Mat
3316: Input Parameter:
3317: . mat - the matrix to transpose
3319: Output Parameters:
3320: . B - the transpose
3322: Notes:
3323: If you pass in PETSC_NULL for B an in-place transpose in mat will be done
3325: Level: intermediate
3327: Concepts: matrices^transposing
3329: .seealso: MatMultTranspose(), MatMultTransposeAdd(), MatIsTranspose()
3330: @*/
3331: PetscErrorCode MatTranspose(Mat mat,Mat *B)
3332: {
3338: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3339: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3340: if (!mat->ops->transpose) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3341: MatPreallocated(mat);
3344: (*mat->ops->transpose)(mat,B);
3346: if (B) {PetscObjectStateIncrease((PetscObject)*B);}
3347: return(0);
3348: }
3352: /*@C
3353: MatIsTranspose - Test whether a matrix is another one's transpose,
3354: or its own, in which case it tests symmetry.
3356: Collective on Mat
3358: Input Parameter:
3359: + A - the matrix to test
3360: - B - the matrix to test against, this can equal the first parameter
3362: Output Parameters:
3363: . flg - the result
3365: Notes:
3366: Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm
3367: has a running time of the order of the number of nonzeros; the parallel
3368: test involves parallel copies of the block-offdiagonal parts of the matrix.
3370: Level: intermediate
3372: Concepts: matrices^transposing, matrix^symmetry
3374: .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian()
3375: @*/
3376: PetscErrorCode MatIsTranspose(Mat A,Mat B,PetscReal tol,PetscTruth *flg)
3377: {
3378: PetscErrorCode ierr,(*f)(Mat,Mat,PetscReal,PetscTruth*),(*g)(Mat,Mat,PetscReal,PetscTruth*);
3384: PetscObjectQueryFunction((PetscObject)A,"MatIsTranspose_C",(void (**)(void))&f);
3385: PetscObjectQueryFunction((PetscObject)B,"MatIsTranspose_C",(void (**)(void))&g);
3386: if (f && g) {
3387: if (f==g) {
3388: (*f)(A,B,tol,flg);
3389: } else {
3390: SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Matrices do not have the same comparator for symmetry test");
3391: }
3392: }
3393: return(0);
3394: }
3398: /*@C
3399: MatPermute - Creates a new matrix with rows and columns permuted from the
3400: original.
3402: Collective on Mat
3404: Input Parameters:
3405: + mat - the matrix to permute
3406: . row - row permutation, each processor supplies only the permutation for its rows
3407: - col - column permutation, each processor needs the entire column permutation, that is
3408: this is the same size as the total number of columns in the matrix
3410: Output Parameters:
3411: . B - the permuted matrix
3413: Level: advanced
3415: Concepts: matrices^permuting
3417: .seealso: MatGetOrdering()
3418: @*/
3419: PetscErrorCode MatPermute(Mat mat,IS row,IS col,Mat *B)
3420: {
3429: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3430: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3431: if (!mat->ops->permute) SETERRQ1(PETSC_ERR_SUP,"MatPermute not available for Mat type %s",mat->type_name);
3432: MatPreallocated(mat);
3434: (*mat->ops->permute)(mat,row,col,B);
3435: PetscObjectStateIncrease((PetscObject)*B);
3436: return(0);
3437: }
3441: /*@C
3442: MatPermuteSparsify - Creates a new matrix with rows and columns permuted from the
3443: original and sparsified to the prescribed tolerance.
3445: Collective on Mat
3447: Input Parameters:
3448: + A - The matrix to permute
3449: . band - The half-bandwidth of the sparsified matrix, or PETSC_DECIDE
3450: . frac - The half-bandwidth as a fraction of the total size, or 0.0
3451: . tol - The drop tolerance
3452: . rowp - The row permutation
3453: - colp - The column permutation
3455: Output Parameter:
3456: . B - The permuted, sparsified matrix
3458: Level: advanced
3460: Note:
3461: The default behavior (band = PETSC_DECIDE and frac = 0.0) is to
3462: restrict the half-bandwidth of the resulting matrix to 5% of the
3463: total matrix size.
3465: .keywords: matrix, permute, sparsify
3467: .seealso: MatGetOrdering(), MatPermute()
3468: @*/
3469: PetscErrorCode MatPermuteSparsify(Mat A, PetscInt band, PetscReal frac, PetscReal tol, IS rowp, IS colp, Mat *B)
3470: {
3471: IS irowp, icolp;
3472: PetscInt *rows, *cols;
3473: PetscInt M, N, locRowStart, locRowEnd;
3474: PetscInt nz, newNz;
3475: const PetscInt *cwork;
3476: PetscInt *cnew;
3477: const PetscScalar *vwork;
3478: PetscScalar *vnew;
3479: PetscInt bw, issize;
3480: PetscInt row, locRow, newRow, col, newCol;
3481: PetscErrorCode ierr;
3488: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
3489: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3490: if (!A->ops->permutesparsify) {
3491: MatGetSize(A, &M, &N);
3492: MatGetOwnershipRange(A, &locRowStart, &locRowEnd);
3493: ISGetSize(rowp, &issize);
3494: if (issize != M) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %D for row permutation, should be %D", issize, M);
3495: ISGetSize(colp, &issize);
3496: if (issize != N) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %D for column permutation, should be %D", issize, N);
3497: ISInvertPermutation(rowp, 0, &irowp);
3498: ISGetIndices(irowp, &rows);
3499: ISInvertPermutation(colp, 0, &icolp);
3500: ISGetIndices(icolp, &cols);
3501: PetscMalloc(N * sizeof(PetscInt), &cnew);
3502: PetscMalloc(N * sizeof(PetscScalar), &vnew);
3504: /* Setup bandwidth to include */
3505: if (band == PETSC_DECIDE) {
3506: if (frac <= 0.0)
3507: bw = (PetscInt) (M * 0.05);
3508: else
3509: bw = (PetscInt) (M * frac);
3510: } else {
3511: if (band <= 0) SETERRQ(PETSC_ERR_ARG_WRONG, "Bandwidth must be a positive integer");
3512: bw = band;
3513: }
3515: /* Put values into new matrix */
3516: MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, B);
3517: for(row = locRowStart, locRow = 0; row < locRowEnd; row++, locRow++) {
3518: MatGetRow(A, row, &nz, &cwork, &vwork);
3519: newRow = rows[locRow]+locRowStart;
3520: for(col = 0, newNz = 0; col < nz; col++) {
3521: newCol = cols[cwork[col]];
3522: if ((newCol >= newRow - bw) && (newCol < newRow + bw) && (PetscAbsScalar(vwork[col]) >= tol)) {
3523: cnew[newNz] = newCol;
3524: vnew[newNz] = vwork[col];
3525: newNz++;
3526: }
3527: }
3528: MatSetValues(*B, 1, &newRow, newNz, cnew, vnew, INSERT_VALUES);
3529: MatRestoreRow(A, row, &nz, &cwork, &vwork);
3530: }
3531: PetscFree(cnew);
3532: PetscFree(vnew);
3533: MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
3534: MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
3535: ISRestoreIndices(irowp, &rows);
3536: ISRestoreIndices(icolp, &cols);
3537: ISDestroy(irowp);
3538: ISDestroy(icolp);
3539: } else {
3540: (*A->ops->permutesparsify)(A, band, frac, tol, rowp, colp, B);
3541: }
3542: PetscObjectStateIncrease((PetscObject)*B);
3543: return(0);
3544: }
3548: /*@
3549: MatEqual - Compares two matrices.
3551: Collective on Mat
3553: Input Parameters:
3554: + A - the first matrix
3555: - B - the second matrix
3557: Output Parameter:
3558: . flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.
3560: Level: intermediate
3562: Concepts: matrices^equality between
3563: @*/
3564: PetscErrorCode MatEqual(Mat A,Mat B,PetscTruth *flg)
3565: {
3573: MatPreallocated(B);
3576: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3577: if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3578: 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);
3579: if (!A->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",A->type_name);
3580: if (!B->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",B->type_name);
3581: if (A->ops->equal != B->ops->equal) SETERRQ2(PETSC_ERR_ARG_INCOMP,"A is type: %s\nB is type: %s",A->type_name,B->type_name);
3582: MatPreallocated(A);
3584: (*A->ops->equal)(A,B,flg);
3585: return(0);
3586: }
3590: /*@
3591: MatDiagonalScale - Scales a matrix on the left and right by diagonal
3592: matrices that are stored as vectors. Either of the two scaling
3593: matrices can be PETSC_NULL.
3595: Collective on Mat
3597: Input Parameters:
3598: + mat - the matrix to be scaled
3599: . l - the left scaling vector (or PETSC_NULL)
3600: - r - the right scaling vector (or PETSC_NULL)
3602: Notes:
3603: MatDiagonalScale() computes A = LAR, where
3604: L = a diagonal matrix (stored as a vector), R = a diagonal matrix (stored as a vector)
3606: Level: intermediate
3608: Concepts: matrices^diagonal scaling
3609: Concepts: diagonal scaling of matrices
3611: .seealso: MatScale()
3612: @*/
3613: PetscErrorCode MatDiagonalScale(Mat mat,Vec l,Vec r)
3614: {
3620: if (!mat->ops->diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3623: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3624: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3625: MatPreallocated(mat);
3628: (*mat->ops->diagonalscale)(mat,l,r);
3630: PetscObjectStateIncrease((PetscObject)mat);
3631: return(0);
3632: }
3636: /*@
3637: MatScale - Scales all elements of a matrix by a given number.
3639: Collective on Mat
3641: Input Parameters:
3642: + mat - the matrix to be scaled
3643: - a - the scaling value
3645: Output Parameter:
3646: . mat - the scaled matrix
3648: Level: intermediate
3650: Concepts: matrices^scaling all entries
3652: .seealso: MatDiagonalScale()
3653: @*/
3654: PetscErrorCode MatScale(Mat mat,PetscScalar a)
3655: {
3661: if (!mat->ops->scale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3662: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3663: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3664: MatPreallocated(mat);
3667: (*mat->ops->scale)(mat,a);
3669: PetscObjectStateIncrease((PetscObject)mat);
3670: return(0);
3671: }
3675: /*@
3676: MatNorm - Calculates various norms of a matrix.
3678: Collective on Mat
3680: Input Parameters:
3681: + mat - the matrix
3682: - type - the type of norm, NORM_1, NORM_FROBENIUS, NORM_INFINITY
3684: Output Parameters:
3685: . nrm - the resulting norm
3687: Level: intermediate
3689: Concepts: matrices^norm
3690: Concepts: norm^of matrix
3691: @*/
3692: PetscErrorCode MatNorm(Mat mat,NormType type,PetscReal *nrm)
3693: {
3701: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3702: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3703: if (!mat->ops->norm) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3704: MatPreallocated(mat);
3706: (*mat->ops->norm)(mat,type,nrm);
3707: return(0);
3708: }
3710: /*
3711: This variable is used to prevent counting of MatAssemblyBegin() that
3712: are called from within a MatAssemblyEnd().
3713: */
3714: static PetscInt MatAssemblyEnd_InUse = 0;
3717: /*@
3718: MatAssemblyBegin - Begins assembling the matrix. This routine should
3719: be called after completing all calls to MatSetValues().
3721: Collective on Mat
3723: Input Parameters:
3724: + mat - the matrix
3725: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
3726:
3727: Notes:
3728: MatSetValues() generally caches the values. The matrix is ready to
3729: use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
3730: Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
3731: in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
3732: using the matrix.
3734: Level: beginner
3736: Concepts: matrices^assembling
3738: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled()
3739: @*/
3740: PetscErrorCode MatAssemblyBegin(Mat mat,MatAssemblyType type)
3741: {
3747: MatPreallocated(mat);
3748: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.\nDid you forget to call MatSetUnfactored()?");
3749: if (mat->assembled) {
3750: mat->was_assembled = PETSC_TRUE;
3751: mat->assembled = PETSC_FALSE;
3752: }
3753: if (!MatAssemblyEnd_InUse) {
3755: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
3757: } else {
3758: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
3759: }
3760: return(0);
3761: }
3765: /*@
3766: MatAssembled - Indicates if a matrix has been assembled and is ready for
3767: use; for example, in matrix-vector product.
3769: Collective on Mat
3771: Input Parameter:
3772: . mat - the matrix
3774: Output Parameter:
3775: . assembled - PETSC_TRUE or PETSC_FALSE
3777: Level: advanced
3779: Concepts: matrices^assembled?
3781: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin()
3782: @*/
3783: PetscErrorCode MatAssembled(Mat mat,PetscTruth *assembled)
3784: {
3789: *assembled = mat->assembled;
3790: return(0);
3791: }
3795: /*
3796: Processes command line options to determine if/how a matrix
3797: is to be viewed. Called by MatAssemblyEnd() and MatLoad().
3798: */
3799: PetscErrorCode MatView_Private(Mat mat)
3800: {
3801: PetscErrorCode ierr;
3802: PetscTruth flg1,flg2,flg3,flg4,flg6,flg7,flg8;
3803: static PetscTruth incall = PETSC_FALSE;
3804: #if defined(PETSC_USE_SOCKET_VIEWER)
3805: PetscTruth flg5;
3806: #endif
3809: if (incall) return(0);
3810: incall = PETSC_TRUE;
3811: PetscOptionsBegin(mat->comm,mat->prefix,"Matrix Options","Mat");
3812: PetscOptionsName("-mat_view_info","Information on matrix size","MatView",&flg1);
3813: PetscOptionsName("-mat_view_info_detailed","Nonzeros in the matrix","MatView",&flg2);
3814: PetscOptionsName("-mat_view","Print matrix to stdout","MatView",&flg3);
3815: PetscOptionsName("-mat_view_matlab","Print matrix to stdout in a format Matlab can read","MatView",&flg4);
3816: #if defined(PETSC_USE_SOCKET_VIEWER)
3817: PetscOptionsName("-mat_view_socket","Send matrix to socket (can be read from matlab)","MatView",&flg5);
3818: #endif
3819: PetscOptionsName("-mat_view_binary","Save matrix to file in binary format","MatView",&flg6);
3820: PetscOptionsName("-mat_view_draw","Draw the matrix nonzero structure","MatView",&flg7);
3821: PetscOptionsEnd();
3823: if (flg1) {
3824: PetscViewer viewer;
3826: PetscViewerASCIIGetStdout(mat->comm,&viewer);
3827: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
3828: MatView(mat,viewer);
3829: PetscViewerPopFormat(viewer);
3830: }
3831: if (flg2) {
3832: PetscViewer viewer;
3834: PetscViewerASCIIGetStdout(mat->comm,&viewer);
3835: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL);
3836: MatView(mat,viewer);
3837: PetscViewerPopFormat(viewer);
3838: }
3839: if (flg3) {
3840: PetscViewer viewer;
3842: PetscViewerASCIIGetStdout(mat->comm,&viewer);
3843: MatView(mat,viewer);
3844: }
3845: if (flg4) {
3846: PetscViewer viewer;
3848: PetscViewerASCIIGetStdout(mat->comm,&viewer);
3849: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
3850: MatView(mat,viewer);
3851: PetscViewerPopFormat(viewer);
3852: }
3853: #if defined(PETSC_USE_SOCKET_VIEWER)
3854: if (flg5) {
3855: MatView(mat,PETSC_VIEWER_SOCKET_(mat->comm));
3856: PetscViewerFlush(PETSC_VIEWER_SOCKET_(mat->comm));
3857: }
3858: #endif
3859: if (flg6) {
3860: MatView(mat,PETSC_VIEWER_BINARY_(mat->comm));
3861: PetscViewerFlush(PETSC_VIEWER_BINARY_(mat->comm));
3862: }
3863: if (flg7) {
3864: PetscOptionsHasName(mat->prefix,"-mat_view_contour",&flg8);
3865: if (flg8) {
3866: PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);
3867: }
3868: MatView(mat,PETSC_VIEWER_DRAW_(mat->comm));
3869: PetscViewerFlush(PETSC_VIEWER_DRAW_(mat->comm));
3870: if (flg8) {
3871: PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));
3872: }
3873: }
3874: incall = PETSC_FALSE;
3875: return(0);
3876: }
3880: /*@
3881: MatAssemblyEnd - Completes assembling the matrix. This routine should
3882: be called after MatAssemblyBegin().
3884: Collective on Mat
3886: Input Parameters:
3887: + mat - the matrix
3888: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
3890: Options Database Keys:
3891: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
3892: . -mat_view_info_detailed - Prints more detailed info
3893: . -mat_view - Prints matrix in ASCII format
3894: . -mat_view_matlab - Prints matrix in Matlab format
3895: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
3896: . -display <name> - Sets display name (default is host)
3897: . -draw_pause <sec> - Sets number of seconds to pause after display
3898: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
3899: . -viewer_socket_machine <machine>
3900: . -viewer_socket_port <port>
3901: . -mat_view_binary - save matrix to file in binary format
3902: - -viewer_binary_filename <name>
3904: Notes:
3905: MatSetValues() generally caches the values. The matrix is ready to
3906: use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
3907: Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
3908: in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
3909: using the matrix.
3911: Level: beginner
3913: .seealso: MatAssemblyBegin(), MatSetValues(), PetscDrawOpenX(), MatView(), MatAssembled(), PetscViewerSocketOpen()
3914: @*/
3915: PetscErrorCode MatAssemblyEnd(Mat mat,MatAssemblyType type)
3916: {
3917: PetscErrorCode ierr;
3918: static PetscInt inassm = 0;
3919: PetscTruth flg;
3925: inassm++;
3926: MatAssemblyEnd_InUse++;
3927: if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
3929: if (mat->ops->assemblyend) {
3930: (*mat->ops->assemblyend)(mat,type);
3931: }
3933: } else {
3934: if (mat->ops->assemblyend) {
3935: (*mat->ops->assemblyend)(mat,type);
3936: }
3937: }
3939: /* Flush assembly is not a true assembly */
3940: if (type != MAT_FLUSH_ASSEMBLY) {
3941: mat->assembled = PETSC_TRUE; mat->num_ass++;
3942: }
3943: mat->insertmode = NOT_SET_VALUES;
3944: MatAssemblyEnd_InUse--;
3945: PetscObjectStateIncrease((PetscObject)mat);
3946: if (!mat->symmetric_eternal) {
3947: mat->symmetric_set = PETSC_FALSE;
3948: mat->hermitian_set = PETSC_FALSE;
3949: mat->structurally_symmetric_set = PETSC_FALSE;
3950: }
3951: if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) {
3952: MatView_Private(mat);
3953: PetscOptionsHasName(mat->prefix,"-mat_is_symmetric",&flg);
3954: if (flg) {
3955: PetscReal tol = 0.0;
3956: PetscOptionsGetReal(mat->prefix,"-mat_is_symmetric",&tol,PETSC_NULL);
3957: MatIsSymmetric(mat,tol,&flg);
3958: if (flg) {
3959: PetscPrintf(mat->comm,"Matrix is symmetric (tolerance %G)\n",tol);
3960: } else {
3961: PetscPrintf(mat->comm,"Matrix is not symmetric (tolerance %G)\n",tol);
3962: }
3963: }
3964: }
3965: inassm--;
3966: return(0);
3967: }
3972: /*@
3973: MatCompress - Tries to store the matrix in as little space as
3974: possible. May fail if memory is already fully used, since it
3975: tries to allocate new space.
3977: Collective on Mat
3979: Input Parameters:
3980: . mat - the matrix
3982: Level: advanced
3984: @*/
3985: PetscErrorCode MatCompress(Mat mat)
3986: {
3992: MatPreallocated(mat);
3993: if (mat->ops->compress) {(*mat->ops->compress)(mat);}
3994: return(0);
3995: }
3999: /*@
4000: MatSetOption - Sets a parameter option for a matrix. Some options
4001: may be specific to certain storage formats. Some options
4002: determine how values will be inserted (or added). Sorted,
4003: row-oriented input will generally assemble the fastest. The default
4004: is row-oriented, nonsorted input.
4006: Collective on Mat
4008: Input Parameters:
4009: + mat - the matrix
4010: - option - the option, one of those listed below (and possibly others),
4011: e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR
4013: Options Describing Matrix Structure:
4014: + MAT_SYMMETRIC - symmetric in terms of both structure and value
4015: . MAT_HERMITIAN - transpose is the complex conjugation
4016: . MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure
4017: . MAT_NOT_SYMMETRIC - not symmetric in value
4018: . MAT_NOT_HERMITIAN - transpose is not the complex conjugation
4019: . MAT_NOT_STRUCTURALLY_SYMMETRIC - not symmetric nonzero structure
4020: . MAT_SYMMETRY_ETERNAL - if you would like the symmetry/Hermitian flag
4021: you set to be kept with all future use of the matrix
4022: including after MatAssemblyBegin/End() which could
4023: potentially change the symmetry structure, i.e. you
4024: KNOW the matrix will ALWAYS have the property you set.
4025: - MAT_NOT_SYMMETRY_ETERNAL - if MatAssemblyBegin/End() is called then the
4026: flags you set will be dropped (in case potentially
4027: the symmetry etc was lost).
4029: Options For Use with MatSetValues():
4030: Insert a logically dense subblock, which can be
4031: + MAT_ROW_ORIENTED - row-oriented (default)
4032: . MAT_COLUMN_ORIENTED - column-oriented
4033: . MAT_ROWS_SORTED - sorted by row
4034: . MAT_ROWS_UNSORTED - not sorted by row (default)
4035: . MAT_COLUMNS_SORTED - sorted by column
4036: - MAT_COLUMNS_UNSORTED - not sorted by column (default)
4038: Not these options reflect the data you pass in with MatSetValues(); it has
4039: nothing to do with how the data is stored internally in the matrix
4040: data structure.
4042: When (re)assembling a matrix, we can restrict the input for
4043: efficiency/debugging purposes. These options include
4044: + MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be
4045: allowed if they generate a new nonzero
4046: . MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed
4047: . MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if
4048: they generate a nonzero in a new diagonal (for block diagonal format only)
4049: . MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only)
4050: . MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries
4051: . MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry
4052: - MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly
4054: Notes:
4055: Some options are relevant only for particular matrix types and
4056: are thus ignored by others. Other options are not supported by
4057: certain matrix types and will generate an error message if set.
4059: If using a Fortran 77 module to compute a matrix, one may need to
4060: use the column-oriented option (or convert to the row-oriented
4061: format).
4063: MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion
4064: that would generate a new entry in the nonzero structure is instead
4065: ignored. Thus, if memory has not alredy been allocated for this particular
4066: data, then the insertion is ignored. For dense matrices, in which
4067: the entire array is allocated, no entries are ever ignored.
4068: Set after the first MatAssemblyEnd()
4070: MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion
4071: that would generate a new entry in the nonzero structure instead produces
4072: an error. (Currently supported for AIJ and BAIJ formats only.)
4073: This is a useful flag when using SAME_NONZERO_PATTERN in calling
4074: KSPSetOperators() to ensure that the nonzero pattern truely does
4075: remain unchanged. Set after the first MatAssemblyEnd()
4077: MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion
4078: that would generate a new entry that has not been preallocated will
4079: instead produce an error. (Currently supported for AIJ and BAIJ formats
4080: only.) This is a useful flag when debugging matrix memory preallocation.
4082: MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for
4083: other processors should be dropped, rather than stashed.
4084: This is useful if you know that the "owning" processor is also
4085: always generating the correct matrix entries, so that PETSc need
4086: not transfer duplicate entries generated on another processor.
4087:
4088: MAT_USE_HASH_TABLE indicates that a hash table be used to improve the
4089: searches during matrix assembly. When this flag is set, the hash table
4090: is created during the first Matrix Assembly. This hash table is
4091: used the next time through, during MatSetVaules()/MatSetVaulesBlocked()
4092: to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag
4093: should be used with MAT_USE_HASH_TABLE flag. This option is currently
4094: supported by MATMPIBAIJ format only.
4096: MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries
4097: are kept in the nonzero structure
4099: MAT_IGNORE_ZERO_ENTRIES - for AIJ/IS matrices this will stop zero values from creating
4100: a zero location in the matrix
4102: MAT_USE_INODES - indicates using inode version of the code - works with AIJ and
4103: ROWBS matrix types
4105: MAT_DO_NOT_USE_INODES - indicates not using inode version of the code - works
4106: with AIJ and ROWBS matrix types (database option "-mat_no_inode")
4108: Level: intermediate
4110: Concepts: matrices^setting options
4112: @*/
4113: PetscErrorCode MatSetOption(Mat mat,MatOption op)
4114: {
4120: MatPreallocated(mat);
4121: switch (op) {
4122: case MAT_SYMMETRIC:
4123: mat->symmetric = PETSC_TRUE;
4124: mat->structurally_symmetric = PETSC_TRUE;
4125: mat->symmetric_set = PETSC_TRUE;
4126: mat->structurally_symmetric_set = PETSC_TRUE;
4127: break;
4128: case MAT_HERMITIAN:
4129: mat->hermitian = PETSC_TRUE;
4130: mat->structurally_symmetric = PETSC_TRUE;
4131: mat->hermitian_set = PETSC_TRUE;
4132: mat->structurally_symmetric_set = PETSC_TRUE;
4133: break;
4134: case MAT_STRUCTURALLY_SYMMETRIC:
4135: mat->structurally_symmetric = PETSC_TRUE;
4136: mat->structurally_symmetric_set = PETSC_TRUE;
4137: break;
4138: case MAT_NOT_SYMMETRIC:
4139: mat->symmetric = PETSC_FALSE;
4140: mat->symmetric_set = PETSC_TRUE;
4141: break;
4142: case MAT_NOT_HERMITIAN:
4143: mat->hermitian = PETSC_FALSE;
4144: mat->hermitian_set = PETSC_TRUE;
4145: break;
4146: case MAT_NOT_STRUCTURALLY_SYMMETRIC:
4147: mat->structurally_symmetric = PETSC_FALSE;
4148: mat->structurally_symmetric_set = PETSC_TRUE;
4149: break;
4150: case MAT_SYMMETRY_ETERNAL:
4151: mat->symmetric_eternal = PETSC_TRUE;
4152: break;
4153: case MAT_NOT_SYMMETRY_ETERNAL:
4154: mat->symmetric_eternal = PETSC_FALSE;
4155: break;
4156: default:
4157: break;
4158: }
4159: if (mat->ops->setoption) {
4160: (*mat->ops->setoption)(mat,op);
4161: }
4162: return(0);
4163: }
4167: /*@
4168: MatZeroEntries - Zeros all entries of a matrix. For sparse matrices
4169: this routine retains the old nonzero structure.
4171: Collective on Mat
4173: Input Parameters:
4174: . mat - the matrix
4176: Level: intermediate
4178: Concepts: matrices^zeroing
4180: .seealso: MatZeroRows()
4181: @*/
4182: PetscErrorCode MatZeroEntries(Mat mat)
4183: {
4189: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4190: if (mat->insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for matrices where you have set values but not yet assembled");
4191: if (!mat->ops->zeroentries) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4192: MatPreallocated(mat);
4195: (*mat->ops->zeroentries)(mat);
4197: PetscObjectStateIncrease((PetscObject)mat);
4198: return(0);
4199: }
4203: /*@C
4204: MatZeroRows - Zeros all entries (except possibly the main diagonal)
4205: of a set of rows of a matrix.
4207: Collective on Mat
4209: Input Parameters:
4210: + mat - the matrix
4211: . numRows - the number of rows to remove
4212: . rows - the global row indices
4213: - diag - value put in all diagonals of eliminated rows
4215: Notes:
4216: For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
4217: but does not release memory. For the dense and block diagonal
4218: formats this does not alter the nonzero structure.
4220: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
4221: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
4222: merely zeroed.
4224: The user can set a value in the diagonal entry (or for the AIJ and
4225: row formats can optionally remove the main diagonal entry from the
4226: nonzero structure as well, by passing 0.0 as the final argument).
4228: For the parallel case, all processes that share the matrix (i.e.,
4229: those in the communicator used for matrix creation) MUST call this
4230: routine, regardless of whether any rows being zeroed are owned by
4231: them.
4233: Each processor should list the rows that IT wants zeroed
4235: Level: intermediate
4237: Concepts: matrices^zeroing rows
4239: .seealso: MatZeroRowsIS(), MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
4240: @*/
4241: PetscErrorCode MatZeroRows(Mat mat,PetscInt numRows,const PetscInt rows[],PetscScalar diag)
4242: {
4249: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4250: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4251: if (!mat->ops->zerorows) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4252: MatPreallocated(mat);
4254: (*mat->ops->zerorows)(mat,numRows,rows,diag);
4255: MatView_Private(mat);
4256: PetscObjectStateIncrease((PetscObject)mat);
4257: return(0);
4258: }
4262: /*@C
4263: MatZeroRowsIS - Zeros all entries (except possibly the main diagonal)
4264: of a set of rows of a matrix.
4266: Collective on Mat
4268: Input Parameters:
4269: + mat - the matrix
4270: . is - index set of rows to remove
4271: - diag - value put in all diagonals of eliminated rows
4273: Notes:
4274: For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
4275: but does not release memory. For the dense and block diagonal
4276: formats this does not alter the nonzero structure.
4278: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
4279: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
4280: merely zeroed.
4282: The user can set a value in the diagonal entry (or for the AIJ and
4283: row formats can optionally remove the main diagonal entry from the
4284: nonzero structure as well, by passing 0.0 as the final argument).
4286: For the parallel case, all processes that share the matrix (i.e.,
4287: those in the communicator used for matrix creation) MUST call this
4288: routine, regardless of whether any rows being zeroed are owned by
4289: them.
4291: Each processor should list the rows that IT wants zeroed
4293: Level: intermediate
4295: Concepts: matrices^zeroing rows
4297: .seealso: MatZeroRows(), MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
4298: @*/
4299: PetscErrorCode MatZeroRowsIS(Mat mat,IS is,PetscScalar diag)
4300: {
4301: PetscInt numRows;
4302: PetscInt *rows;
4309: ISGetLocalSize(is,&numRows);
4310: ISGetIndices(is,&rows);
4311: MatZeroRows(mat,numRows,rows,diag);
4312: ISRestoreIndices(is,&rows);
4313: return(0);
4314: }
4318: /*@C
4319: MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
4320: of a set of rows of a matrix; using local numbering of rows.
4322: Collective on Mat
4324: Input Parameters:
4325: + mat - the matrix
4326: . numRows - the number of rows to remove
4327: . rows - the global row indices
4328: - diag - value put in all diagonals of eliminated rows
4330: Notes:
4331: Before calling MatZeroRowsLocal(), the user must first set the
4332: local-to-global mapping by calling MatSetLocalToGlobalMapping().
4334: For the AIJ matrix formats this removes the old nonzero structure,
4335: but does not release memory. For the dense and block diagonal
4336: formats this does not alter the nonzero structure.
4338: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
4339: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
4340: merely zeroed.
4342: The user can set a value in the diagonal entry (or for the AIJ and
4343: row formats can optionally remove the main diagonal entry from the
4344: nonzero structure as well, by passing 0.0 as the final argument).
4346: Level: intermediate
4348: Concepts: matrices^zeroing
4350: .seealso: MatZeroRows(), MatZeroRowsLocalIS(), MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
4351: @*/
4352: PetscErrorCode MatZeroRowsLocal(Mat mat,PetscInt numRows,const PetscInt rows[],PetscScalar diag)
4353: {
4360: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4361: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4362: MatPreallocated(mat);
4364: if (mat->ops->zerorowslocal) {
4365: (*mat->ops->zerorowslocal)(mat,numRows,rows,diag);
4366: } else {
4367: IS is, newis;
4368: PetscInt *newRows;
4370: if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Need to provide local to global mapping to matrix first");
4371: ISCreateGeneral(PETSC_COMM_SELF,numRows,rows,&is);
4372: ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);
4373: ISGetIndices(newis,&newRows);
4374: (*mat->ops->zerorows)(mat,numRows,newRows,diag);
4375: ISRestoreIndices(newis,&newRows);
4376: ISDestroy(newis);
4377: ISDestroy(is);
4378: }
4379: PetscObjectStateIncrease((PetscObject)mat);
4380: return(0);
4381: }
4385: /*@C
4386: MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
4387: of a set of rows of a matrix; using local numbering of rows.
4389: Collective on Mat
4391: Input Parameters:
4392: + mat - the matrix
4393: . is - index set of rows to remove
4394: - diag - value put in all diagonals of eliminated rows
4396: Notes:
4397: Before calling MatZeroRowsLocal(), the user must first set the
4398: local-to-global mapping by calling MatSetLocalToGlobalMapping().
4400: For the AIJ matrix formats this removes the old nonzero structure,
4401: but does not release memory. For the dense and block diagonal
4402: formats this does not alter the nonzero structure.
4404: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
4405: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
4406: merely zeroed.
4408: The user can set a value in the diagonal entry (or for the AIJ and
4409: row formats can optionally remove the main diagonal entry from the
4410: nonzero structure as well, by passing 0.0 as the final argument).
4412: Level: intermediate
4414: Concepts: matrices^zeroing
4416: .seealso: MatZeroRows(), MatZeroRowsLocal(), MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
4417: @*/
4418: PetscErrorCode MatZeroRowsLocalIS(Mat mat,IS is,PetscScalar diag)
4419: {
4421: PetscInt numRows;
4422: PetscInt *rows;
4428: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4429: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4430: MatPreallocated(mat);
4432: ISGetLocalSize(is,&numRows);
4433: ISGetIndices(is,&rows);
4434: MatZeroRowsLocal(mat,numRows,rows,diag);
4435: ISRestoreIndices(is,&rows);
4436: return(0);
4437: }
4441: /*@
4442: MatGetSize - Returns the numbers of rows and columns in a matrix.
4444: Not Collective
4446: Input Parameter:
4447: . mat - the matrix
4449: Output Parameters:
4450: + m - the number of global rows
4451: - n - the number of global columns
4453: Note: both output parameters can be PETSC_NULL on input.
4455: Level: beginner
4457: Concepts: matrices^size
4459: .seealso: MatGetLocalSize()
4460: @*/
4461: PetscErrorCode MatGetSize(Mat mat,PetscInt *m,PetscInt* n)
4462: {
4465: if (m) *m = mat->rmap.N;
4466: if (n) *n = mat->cmap.N;
4467: return(0);
4468: }
4472: /*@
4473: MatGetLocalSize - Returns the number of rows and columns in a matrix
4474: stored locally. This information may be implementation dependent, so
4475: use with care.
4477: Not Collective
4479: Input Parameters:
4480: . mat - the matrix
4482: Output Parameters:
4483: + m - the number of local rows
4484: - n - the number of local columns
4486: Note: both output parameters can be PETSC_NULL on input.
4488: Level: beginner
4490: Concepts: matrices^local size
4492: .seealso: MatGetSize()
4493: @*/
4494: PetscErrorCode MatGetLocalSize(Mat mat,PetscInt *m,PetscInt* n)
4495: {
4500: if (m) *m = mat->rmap.n;
4501: if (n) *n = mat->cmap.n;
4502: return(0);
4503: }
4508: /*@
4509: MatGetOwnershipRange - Returns the range of matrix rows owned by
4510: this processor, assuming that the matrix is laid out with the first
4511: n1 rows on the first processor, the next n2 rows on the second, etc.
4512: For certain parallel layouts this range may not be well defined.
4514: Not Collective
4516: Input Parameters:
4517: . mat - the matrix
4519: Output Parameters:
4520: + m - the global index of the first local row
4521: - n - one more than the global index of the last local row
4523: Note: both output parameters can be PETSC_NULL on input.
4525: Level: beginner
4527: Concepts: matrices^row ownership
4529: .seealso: MatGetOwnershipRanges()
4531: @*/
4532: PetscErrorCode MatGetOwnershipRange(Mat mat,PetscInt *m,PetscInt* n)
4533: {
4541: MatPreallocated(mat);
4542: if (m) *m = mat->rmap.rstart;
4543: if (n) *n = mat->rmap.rend;
4544: return(0);
4545: }
4549: /*@C
4550: MatGetOwnershipRanges - Returns the range of matrix rows owned by
4551: each process
4553: Not Collective
4555: Input Parameters:
4556: . mat - the matrix
4558: Output Parameters:
4559: . ranges - start of each processors portion plus one more then the total length at the end
4561: Level: beginner
4563: Concepts: matrices^row ownership
4565: .seealso: MatGetOwnershipRange()
4567: @*/
4568: PetscErrorCode MatGetOwnershipRanges(Mat mat,const PetscInt **ranges)
4569: {
4575: PetscMapGetGlobalRange(&mat->rmap,ranges);
4576: return(0);
4577: }
4581: /*@
4582: MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
4583: Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric()
4584: to complete the factorization.
4586: Collective on Mat
4588: Input Parameters:
4589: + mat - the matrix
4590: . row - row permutation
4591: . column - column permutation
4592: - info - structure containing
4593: $ levels - number of levels of fill.
4594: $ expected fill - as ratio of original fill.
4595: $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices
4596: missing diagonal entries)
4598: Output Parameters:
4599: . fact - new matrix that has been symbolically factored
4601: Notes:
4602: See the users manual for additional information about
4603: choosing the fill factor for better efficiency.
4605: Most users should employ the simplified KSP interface for linear solvers
4606: instead of working directly with matrix algebra routines such as this.
4607: See, e.g., KSPCreate().
4609: Level: developer
4611: Concepts: matrices^symbolic LU factorization
4612: Concepts: matrices^factorization
4613: Concepts: LU^symbolic factorization
4615: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
4616: MatGetOrdering(), MatFactorInfo
4618: @*/
4619: PetscErrorCode MatILUFactorSymbolic(Mat mat,IS row,IS col,MatFactorInfo *info,Mat *fact)
4620: {
4630: if (info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels of fill negative %D",(PetscInt)info->levels);
4631: if (info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %G",info->fill);
4632: if (!mat->ops->ilufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic ILU",mat->type_name);
4633: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4634: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4635: MatPreallocated(mat);
4638: (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);
4640: return(0);
4641: }
4645: /*@
4646: MatICCFactorSymbolic - Performs symbolic incomplete
4647: Cholesky factorization for a symmetric matrix. Use
4648: MatCholeskyFactorNumeric() to complete the factorization.
4650: Collective on Mat
4652: Input Parameters:
4653: + mat - the matrix
4654: . perm - row and column permutation
4655: - info - structure containing
4656: $ levels - number of levels of fill.
4657: $ expected fill - as ratio of original fill.
4659: Output Parameter:
4660: . fact - the factored matrix
4662: Notes:
4663: Most users should employ the KSP interface for linear solvers
4664: instead of working directly with matrix algebra routines such as this.
4665: See, e.g., KSPCreate().
4667: Level: developer
4669: Concepts: matrices^symbolic incomplete Cholesky factorization
4670: Concepts: matrices^factorization
4671: Concepts: Cholsky^symbolic factorization
4673: .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
4674: @*/
4675: PetscErrorCode MatICCFactorSymbolic(Mat mat,IS perm,MatFactorInfo *info,Mat *fact)
4676: {
4685: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4686: if (info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels negative %D",(PetscInt) info->levels);
4687: if (info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %G",info->fill);
4688: if (!mat->ops->iccfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic ICC",mat->type_name);
4689: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4690: MatPreallocated(mat);
4693: (*mat->ops->iccfactorsymbolic)(mat,perm,info,fact);
4695: return(0);
4696: }
4700: /*@C
4701: MatGetArray - Returns a pointer to the element values in the matrix.
4702: The result of this routine is dependent on the underlying matrix data
4703: structure, and may not even work for certain matrix types. You MUST
4704: call MatRestoreArray() when you no longer need to access the array.
4706: Not Collective
4708: Input Parameter:
4709: . mat - the matrix
4711: Output Parameter:
4712: . v - the location of the values
4715: Fortran Note:
4716: This routine is used differently from Fortran, e.g.,
4717: .vb
4718: Mat mat
4719: PetscScalar mat_array(1)
4720: PetscOffset i_mat
4721: PetscErrorCode ierr
4722: call MatGetArray(mat,mat_array,i_mat,ierr)
4724: C Access first local entry in matrix; note that array is
4725: C treated as one dimensional
4726: value = mat_array(i_mat + 1)
4728: [... other code ...]
4729: call MatRestoreArray(mat,mat_array,i_mat,ierr)
4730: .ve
4732: See the Fortran chapter of the users manual and
4733: petsc/src/mat/examples/tests for details.
4735: Level: advanced
4737: Concepts: matrices^access array
4739: .seealso: MatRestoreArray(), MatGetArrayF90()
4740: @*/
4741: PetscErrorCode MatGetArray(Mat mat,PetscScalar *v[])
4742: {
4749: if (!mat->ops->getarray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4750: MatPreallocated(mat);
4751: (*mat->ops->getarray)(mat,v);
4752: CHKMEMQ;
4753: return(0);
4754: }
4758: /*@C
4759: MatRestoreArray - Restores the matrix after MatGetArray() has been called.
4761: Not Collective
4763: Input Parameter:
4764: + mat - the matrix
4765: - v - the location of the values
4767: Fortran Note:
4768: This routine is used differently from Fortran, e.g.,
4769: .vb
4770: Mat mat
4771: PetscScalar mat_array(1)
4772: PetscOffset i_mat
4773: PetscErrorCode ierr
4774: call MatGetArray(mat,mat_array,i_mat,ierr)
4776: C Access first local entry in matrix; note that array is
4777: C treated as one dimensional
4778: value = mat_array(i_mat + 1)
4780: [... other code ...]
4781: call MatRestoreArray(mat,mat_array,i_mat,ierr)
4782: .ve
4784: See the Fortran chapter of the users manual and
4785: petsc/src/mat/examples/tests for details
4787: Level: advanced
4789: .seealso: MatGetArray(), MatRestoreArrayF90()
4790: @*/
4791: PetscErrorCode MatRestoreArray(Mat mat,PetscScalar *v[])
4792: {
4799: #if defined(PETSC_USE_DEBUG)
4800: CHKMEMQ;
4801: #endif
4802: if (!mat->ops->restorearray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4803: (*mat->ops->restorearray)(mat,v);
4804: PetscObjectStateIncrease((PetscObject)mat);
4805: return(0);
4806: }
4810: /*@C
4811: MatGetSubMatrices - Extracts several submatrices from a matrix. If submat
4812: points to an array of valid matrices, they may be reused to store the new
4813: submatrices.
4815: Collective on Mat
4817: Input Parameters:
4818: + mat - the matrix
4819: . n - the number of submatrixes to be extracted (on this processor, may be zero)
4820: . irow, icol - index sets of rows and columns to extract
4821: - scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
4823: Output Parameter:
4824: . submat - the array of submatrices
4826: Notes:
4827: MatGetSubMatrices() can extract only sequential submatrices
4828: (from both sequential and parallel matrices). Use MatGetSubMatrix()
4829: to extract a parallel submatrix.
4831: When extracting submatrices from a parallel matrix, each processor can
4832: form a different submatrix by setting the rows and columns of its
4833: individual index sets according to the local submatrix desired.
4835: When finished using the submatrices, the user should destroy
4836: them with MatDestroyMatrices().
4838: MAT_REUSE_MATRIX can only be used when the nonzero structure of the
4839: original matrix has not changed from that last call to MatGetSubMatrices().
4841: This routine creates the matrices in submat; you should NOT create them before
4842: calling it. It also allocates the array of matrix pointers submat.
4844: For BAIJ matrices the index sets must respect the block structure, that is if they
4845: request one row/column in a block, they must request all rows/columns that are in
4846: that block. For example, if the block size is 2 you cannot request just row 0 and
4847: column 0.
4849: Fortran Note:
4850: The Fortran interface is slightly different from that given below; it
4851: requires one to pass in as submat a Mat (integer) array of size at least m.
4853: Level: advanced
4855: Concepts: matrices^accessing submatrices
4856: Concepts: submatrices
4858: .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal()
4859: @*/
4860: PetscErrorCode MatGetSubMatrices(Mat mat,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *submat[])
4861: {
4863: PetscInt i;
4864: PetscTruth eq;
4869: if (n) {
4874: }
4876: if (n && scall == MAT_REUSE_MATRIX) {
4879: }
4880: if (!mat->ops->getsubmatrices) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4881: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4882: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4883: MatPreallocated(mat);
4886: (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);
4888: for (i=0; i<n; i++) {
4889: if (mat->symmetric || mat->structurally_symmetric || mat->hermitian) {
4890: ISEqual(irow[i],icol[i],&eq);
4891: if (eq) {
4892: if (mat->symmetric){
4893: MatSetOption((*submat)[i],MAT_SYMMETRIC);
4894: } else if (mat->hermitian) {
4895: MatSetOption((*submat)[i],MAT_HERMITIAN);
4896: } else if (mat->structurally_symmetric) {
4897: MatSetOption((*submat)[i],MAT_STRUCTURALLY_SYMMETRIC);
4898: }
4899: }
4900: }
4901: }
4902: return(0);
4903: }
4907: /*@C
4908: MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices().
4910: Collective on Mat
4912: Input Parameters:
4913: + n - the number of local matrices
4914: - mat - the matrices (note that this is a pointer to the array of matrices, just to match the calling
4915: sequence of MatGetSubMatrices())
4917: Level: advanced
4919: Notes: Frees not only the matrices, but also the array that contains the matrices
4921: .seealso: MatGetSubMatrices()
4922: @*/
4923: PetscErrorCode MatDestroyMatrices(PetscInt n,Mat *mat[])
4924: {
4926: PetscInt i;
4929: if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of matrices %D",n);
4931: for (i=0; i<n; i++) {
4932: MatDestroy((*mat)[i]);
4933: }
4934: /* memory is allocated even if n = 0 */
4935: PetscFree(*mat);
4936: return(0);
4937: }
4941: /*@
4942: MatIncreaseOverlap - Given a set of submatrices indicated by index sets,
4943: replaces the index sets by larger ones that represent submatrices with
4944: additional overlap.
4946: Collective on Mat
4948: Input Parameters:
4949: + mat - the matrix
4950: . n - the number of index sets
4951: . is - the array of index sets (these index sets will changed during the call)
4952: - ov - the additional overlap requested
4954: Level: developer
4956: Concepts: overlap
4957: Concepts: ASM^computing overlap
4959: .seealso: MatGetSubMatrices()
4960: @*/
4961: PetscErrorCode MatIncreaseOverlap(Mat mat,PetscInt n,IS is[],PetscInt ov)
4962: {
4968: if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more domains, you have %D",n);
4969: if (n) {
4972: }
4973: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4974: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4975: MatPreallocated(mat);
4977: if (!ov) return(0);
4978: if (!mat->ops->increaseoverlap) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4980: (*mat->ops->increaseoverlap)(mat,n,is,ov);
4982: return(0);
4983: }
4987: /*@
4988: MatGetBlockSize - Returns the matrix block size; useful especially for the
4989: block row and block diagonal formats.
4990:
4991: Not Collective
4993: Input Parameter:
4994: . mat - the matrix
4996: Output Parameter:
4997: . bs - block size
4999: Notes:
5000: Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG.
5001: Block row formats are MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ
5003: Level: intermediate
5005: Concepts: matrices^block size
5007: .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag()
5008: @*/
5009: PetscErrorCode MatGetBlockSize(Mat mat,PetscInt *bs)
5010: {
5017: MatPreallocated(mat);
5018: *bs = mat->rmap.bs;
5019: return(0);
5020: }
5024: /*@
5025: MatSetBlockSize - Sets the matrix block size; for many matrix types you
5026: cannot use this and MUST set the blocksize when you preallocate the matrix
5027:
5028: Not Collective
5030: Input Parameters:
5031: + mat - the matrix
5032: - bs - block size
5034: Notes:
5035: Only works for shell and AIJ matrices
5037: Level: intermediate
5039: Concepts: matrices^block size
5041: .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag(), MatGetBlockSize()
5042: @*/
5043: PetscErrorCode MatSetBlockSize(Mat mat,PetscInt bs)
5044: {
5050: MatPreallocated(mat);
5051: if (mat->ops->setblocksize) {
5052: mat->rmap.bs = bs;
5053: (*mat->ops->setblocksize)(mat,bs);
5054: } else {
5055: SETERRQ1(PETSC_ERR_ARG_INCOMP,"Cannot set the blocksize for matrix type %s",mat->type_name);
5056: }
5057: return(0);
5058: }
5062: /*@C
5063: MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices.
5065: Collective on Mat
5067: Input Parameters:
5068: + mat - the matrix
5069: . shift - 0 or 1 indicating we want the indices starting at 0 or 1
5070: - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
5071: symmetrized
5073: Output Parameters:
5074: + n - number of rows in the (possibly compressed) matrix
5075: . ia - the row pointers
5076: . ja - the column indices
5077: - done - indicates if the routine actually worked and returned appropriate ia[] and ja[] arrays; callers
5078: are responsible for handling the case when done == PETSC_FALSE and ia and ja are not set
5080: Level: developer
5082: Notes: You CANNOT change any of the ia[] or ja[] values.
5084: Use MatRestoreRowIJ() when you are finished accessing the ia[] and ja[] values
5086: .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
5087: @*/
5088: PetscErrorCode MatGetRowIJ(Mat mat,PetscInt shift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt* ja[],PetscTruth *done)
5089: {
5099: MatPreallocated(mat);
5100: if (!mat->ops->getrowij) *done = PETSC_FALSE;
5101: else {
5102: *done = PETSC_TRUE;
5103: (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);
5104: }
5105: return(0);
5106: }
5110: /*@C
5111: MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices.
5113: Collective on Mat
5115: Input Parameters:
5116: + mat - the matrix
5117: . shift - 1 or zero indicating we want the indices starting at 0 or 1
5118: - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
5119: symmetrized
5121: Output Parameters:
5122: + n - number of columns in the (possibly compressed) matrix
5123: . ia - the column pointers
5124: . ja - the row indices
5125: - done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
5127: Level: developer
5129: .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
5130: @*/
5131: PetscErrorCode MatGetColumnIJ(Mat mat,PetscInt shift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt* ja[],PetscTruth *done)
5132: {
5142: MatPreallocated(mat);
5143: if (!mat->ops->getcolumnij) *done = PETSC_FALSE;
5144: else {
5145: *done = PETSC_TRUE;
5146: (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);
5147: }
5148: return(0);
5149: }
5153: /*@C
5154: MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with
5155: MatGetRowIJ().
5157: Collective on Mat
5159: Input Parameters:
5160: + mat - the matrix
5161: . shift - 1 or zero indicating we want the indices starting at 0 or 1
5162: - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
5163: symmetrized
5165: Output Parameters:
5166: + n - size of (possibly compressed) matrix
5167: . ia - the row pointers
5168: . ja - the column indices
5169: - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
5171: Level: developer
5173: .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
5174: @*/
5175: PetscErrorCode MatRestoreRowIJ(Mat mat,PetscInt shift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt* ja[],PetscTruth *done)
5176: {
5185: MatPreallocated(mat);
5187: if (!mat->ops->restorerowij) *done = PETSC_FALSE;
5188: else {
5189: *done = PETSC_TRUE;
5190: (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);
5191: }
5192: return(0);
5193: }
5197: /*@C
5198: MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with
5199: MatGetColumnIJ().
5201: Collective on Mat
5203: Input Parameters:
5204: + mat - the matrix
5205: . shift - 1 or zero indicating we want the indices starting at 0 or 1
5206: - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
5207: symmetrized
5209: Output Parameters:
5210: + n - size of (possibly compressed) matrix
5211: . ia - the column pointers
5212: . ja - the row indices
5213: - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
5215: Level: developer
5217: .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
5218: @*/
5219: PetscErrorCode MatRestoreColumnIJ(Mat mat,PetscInt shift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt* ja[],PetscTruth *done)
5220: {
5229: MatPreallocated(mat);
5231: if (!mat->ops->restorecolumnij) *done = PETSC_FALSE;
5232: else {
5233: *done = PETSC_TRUE;
5234: (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);
5235: }
5236: return(0);
5237: }
5241: /*@C
5242: MatColoringPatch -Used inside matrix coloring routines that
5243: use MatGetRowIJ() and/or MatGetColumnIJ().
5245: Collective on Mat
5247: Input Parameters:
5248: + mat - the matrix
5249: . ncolors - max color value
5250: . n - number of entries in colorarray
5251: - colorarray - array indicating color for each column
5253: Output Parameters:
5254: . iscoloring - coloring generated using colorarray information
5256: Level: developer
5258: .seealso: MatGetRowIJ(), MatGetColumnIJ()
5260: @*/
5261: PetscErrorCode MatColoringPatch(Mat mat,PetscInt ncolors,PetscInt n,ISColoringValue colorarray[],ISColoring *iscoloring)
5262: {
5270: MatPreallocated(mat);
5272: if (!mat->ops->coloringpatch){
5273: ISColoringCreate(mat->comm,ncolors,n,colorarray,iscoloring);
5274: } else {
5275: (*mat->ops->coloringpatch)(mat,ncolors,n,colorarray,iscoloring);
5276: }
5277: return(0);
5278: }
5283: /*@
5284: MatSetUnfactored - Resets a factored matrix to be treated as unfactored.
5286: Collective on Mat
5288: Input Parameter:
5289: . mat - the factored matrix to be reset
5291: Notes:
5292: This routine should be used only with factored matrices formed by in-place
5293: factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE
5294: format). This option can save memory, for example, when solving nonlinear
5295: systems with a matrix-free Newton-Krylov method and a matrix-based, in-place
5296: ILU(0) preconditioner.
5298: Note that one can specify in-place ILU(0) factorization by calling
5299: .vb
5300: PCType(pc,PCILU);
5301: PCFactorSeUseInPlace(pc);
5302: .ve
5303: or by using the options -pc_type ilu -pc_factor_in_place
5305: In-place factorization ILU(0) can also be used as a local
5306: solver for the blocks within the block Jacobi or additive Schwarz
5307: methods (runtime option: -sub_pc_factor_in_place). See the discussion
5308: of these preconditioners in the users manual for details on setting
5309: local solver options.
5311: Most users should employ the simplified KSP interface for linear solvers
5312: instead of working directly with matrix algebra routines such as this.
5313: See, e.g., KSPCreate().
5315: Level: developer
5317: .seealso: PCFactorSetUseInPlace()
5319: Concepts: matrices^unfactored
5321: @*/
5322: PetscErrorCode MatSetUnfactored(Mat mat)
5323: {
5329: MatPreallocated(mat);
5330: mat->factor = 0;
5331: if (!mat->ops->setunfactored) return(0);
5332: (*mat->ops->setunfactored)(mat);
5333: return(0);
5334: }
5336: /*MC
5337: MatGetArrayF90 - Accesses a matrix array from Fortran90.
5339: Synopsis:
5340: MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
5342: Not collective
5344: Input Parameter:
5345: . x - matrix
5347: Output Parameters:
5348: + xx_v - the Fortran90 pointer to the array
5349: - ierr - error code
5351: Example of Usage:
5352: .vb
5353: PetscScalar, pointer xx_v(:)
5354: ....
5355: call MatGetArrayF90(x,xx_v,ierr)
5356: a = xx_v(3)
5357: call MatRestoreArrayF90(x,xx_v,ierr)
5358: .ve
5360: Notes:
5361: Not yet supported for all F90 compilers
5363: Level: advanced
5365: .seealso: MatRestoreArrayF90(), MatGetArray(), MatRestoreArray()
5367: Concepts: matrices^accessing array
5369: M*/
5371: /*MC
5372: MatRestoreArrayF90 - Restores a matrix array that has been
5373: accessed with MatGetArrayF90().
5375: Synopsis:
5376: MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
5378: Not collective
5380: Input Parameters:
5381: + x - matrix
5382: - xx_v - the Fortran90 pointer to the array
5384: Output Parameter:
5385: . ierr - error code
5387: Example of Usage:
5388: .vb
5389: PetscScalar, pointer xx_v(:)
5390: ....
5391: call MatGetArrayF90(x,xx_v,ierr)
5392: a = xx_v(3)
5393: call MatRestoreArrayF90(x,xx_v,ierr)
5394: .ve
5395:
5396: Notes:
5397: Not yet supported for all F90 compilers
5399: Level: advanced
5401: .seealso: MatGetArrayF90(), MatGetArray(), MatRestoreArray()
5403: M*/
5408: /*@
5409: MatGetSubMatrix - Gets a single submatrix on the same number of processors
5410: as the original matrix.
5412: Collective on Mat
5414: Input Parameters:
5415: + mat - the original matrix
5416: . isrow - rows this processor should obtain
5417: . iscol - columns for all processors you wish to keep
5418: . csize - number of columns "local" to this processor (does nothing for sequential
5419: matrices). This should match the result from VecGetLocalSize(x,...) if you
5420: plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE
5421: - cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
5423: Output Parameter:
5424: . newmat - the new submatrix, of the same type as the old
5426: Level: advanced
5428: Notes: the iscol argument MUST be the same on each processor. You might be
5429: able to create the iscol argument with ISAllGather().
5431: The first time this is called you should use a cll of MAT_INITIAL_MATRIX,
5432: the MatGetSubMatrix() routine will create the newmat for you. Any additional calls
5433: to this routine with a mat of the same nonzero structure and with a call of MAT_REUSE_MATRIX
5434: will reuse the matrix generated the first time. You should call MatDestroy() on newmat when
5435: you are finished using it.
5437: Concepts: matrices^submatrices
5439: .seealso: MatGetSubMatrices(), ISAllGather()
5440: @*/
5441: PetscErrorCode MatGetSubMatrix(Mat mat,IS isrow,IS iscol,PetscInt csize,MatReuse cll,Mat *newmat)
5442: {
5444: PetscMPIInt size;
5445: Mat *local;
5454: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
5455: MatPreallocated(mat);
5456: MPI_Comm_size(mat->comm,&size);
5458: /* if original matrix is on just one processor then use submatrix generated */
5459: if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) {
5460: MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);
5461: return(0);
5462: } else if (!mat->ops->getsubmatrix && size == 1) {
5463: MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);
5464: *newmat = *local;
5465: PetscFree(local);
5466: return(0);
5467: }
5469: if (!mat->ops->getsubmatrix) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
5470: (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);
5471: PetscObjectStateIncrease((PetscObject)*newmat);
5472: return(0);
5473: }
5477: /*@
5478: MatGetSubMatrixRaw - Gets a single submatrix on the same number of processors
5479: as the original matrix.
5481: Collective on Mat
5483: Input Parameters:
5484: + mat - the original matrix
5485: . nrows - the number of rows this processor should obtain
5486: . rows - rows this processor should obtain
5487: . ncols - the number of columns for all processors you wish to keep
5488: . cols - columns for all processors you wish to keep
5489: . csize - number of columns "local" to this processor (does nothing for sequential
5490: matrices). This should match the result from VecGetLocalSize(x,...) if you
5491: plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE
5492: - cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
5494: Output Parameter:
5495: . newmat - the new submatrix, of the same type as the old
5497: Level: advanced
5499: Notes: the iscol argument MUST be the same on each processor. You might be
5500: able to create the iscol argument with ISAllGather().
5502: The first time this is called you should use a cll of MAT_INITIAL_MATRIX,
5503: the MatGetSubMatrix() routine will create the newmat for you. Any additional calls
5504: to this routine with a mat of the same nonzero structure and with a cll of MAT_REUSE_MATRIX
5505: will reuse the matrix generated the first time.
5507: Concepts: matrices^submatrices
5509: .seealso: MatGetSubMatrices(), ISAllGather()
5510: @*/
5511: PetscErrorCode MatGetSubMatrixRaw(Mat mat,PetscInt nrows,const PetscInt rows[],PetscInt ncols,const PetscInt cols[],PetscInt csize,MatReuse cll,Mat *newmat)
5512: {
5513: IS isrow, iscol;
5523: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
5524: MatPreallocated(mat);
5525: ISCreateGeneralWithArray(PETSC_COMM_SELF, nrows, (PetscInt *) rows, &isrow);
5526: ISCreateGeneralWithArray(PETSC_COMM_SELF, ncols, (PetscInt *) cols, &iscol);
5527: MatGetSubMatrix(mat, isrow, iscol, csize, cll, newmat);
5528: ISDestroy(isrow);
5529: ISDestroy(iscol);
5530: return(0);
5531: }
5535: /*@
5536: MatStashSetInitialSize - sets the sizes of the matrix stash, that is
5537: used during the assembly process to store values that belong to
5538: other processors.
5540: Not Collective
5542: Input Parameters:
5543: + mat - the matrix
5544: . size - the initial size of the stash.
5545: - bsize - the initial size of the block-stash(if used).
5547: Options Database Keys:
5548: + -matstash_initial_size <size> or <size0,size1,...sizep-1>
5549: - -matstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>
5551: Level: intermediate
5553: Notes:
5554: The block-stash is used for values set with MatSetValuesBlocked() while
5555: the stash is used for values set with MatSetValues()
5557: Run with the option -info and look for output of the form
5558: MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
5559: to determine the appropriate value, MM, to use for size and
5560: MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
5561: to determine the value, BMM to use for bsize
5563: Concepts: stash^setting matrix size
5564: Concepts: matrices^stash
5566: @*/
5567: PetscErrorCode MatStashSetInitialSize(Mat mat,PetscInt size, PetscInt bsize)
5568: {
5574: MatStashSetInitialSize_Private(&mat->stash,size);
5575: MatStashSetInitialSize_Private(&mat->bstash,bsize);
5576: return(0);
5577: }
5581: /*@
5582: MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of
5583: the matrix
5585: Collective on Mat
5587: Input Parameters:
5588: + mat - the matrix
5589: . x,y - the vectors
5590: - w - where the result is stored
5592: Level: intermediate
5594: Notes:
5595: w may be the same vector as y.
5597: This allows one to use either the restriction or interpolation (its transpose)
5598: matrix to do the interpolation
5600: Concepts: interpolation
5602: .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
5604: @*/
5605: PetscErrorCode MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w)
5606: {
5608: PetscInt M,N;
5616: MatPreallocated(A);
5617: MatGetSize(A,&M,&N);
5618: if (N > M) {
5619: MatMultTransposeAdd(A,x,y,w);
5620: } else {
5621: MatMultAdd(A,x,y,w);
5622: }
5623: return(0);
5624: }
5628: /*@
5629: MatInterpolate - y = A*x or A'*x depending on the shape of
5630: the matrix
5632: Collective on Mat
5634: Input Parameters:
5635: + mat - the matrix
5636: - x,y - the vectors
5638: Level: intermediate
5640: Notes:
5641: This allows one to use either the restriction or interpolation (its transpose)
5642: matrix to do the interpolation
5644: Concepts: matrices^interpolation
5646: .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
5648: @*/
5649: PetscErrorCode MatInterpolate(Mat A,Vec x,Vec y)
5650: {
5652: PetscInt M,N;
5659: MatPreallocated(A);
5660: MatGetSize(A,&M,&N);
5661: if (N > M) {
5662: MatMultTranspose(A,x,y);
5663: } else {
5664: MatMult(A,x,y);
5665: }
5666: return(0);
5667: }
5671: /*@
5672: MatRestrict - y = A*x or A'*x
5674: Collective on Mat
5676: Input Parameters:
5677: + mat - the matrix
5678: - x,y - the vectors
5680: Level: intermediate
5682: Notes:
5683: This allows one to use either the restriction or interpolation (its transpose)
5684: matrix to do the restriction
5686: Concepts: matrices^restriction
5688: .seealso: MatMultAdd(), MatMultTransposeAdd(), MatInterpolate()
5690: @*/
5691: PetscErrorCode MatRestrict(Mat A,Vec x,Vec y)
5692: {
5694: PetscInt M,N;
5701: MatPreallocated(A);
5703: MatGetSize(A,&M,&N);
5704: if (N > M) {
5705: MatMult(A,x,y);
5706: } else {
5707: MatMultTranspose(A,x,y);
5708: }
5709: return(0);
5710: }
5714: /*@C
5715: MatNullSpaceAttach - attaches a null space to a matrix.
5716: This null space will be removed from the resulting vector whenever
5717: MatMult() is called
5719: Collective on Mat
5721: Input Parameters:
5722: + mat - the matrix
5723: - nullsp - the null space object
5725: Level: developer
5727: Notes:
5728: Overwrites any previous null space that may have been attached
5730: Concepts: null space^attaching to matrix
5732: .seealso: MatCreate(), MatNullSpaceCreate()
5733: @*/
5734: PetscErrorCode MatNullSpaceAttach(Mat mat,MatNullSpace nullsp)
5735: {
5742: MatPreallocated(mat);
5744: PetscObjectReference((PetscObject)nullsp);
5745: if (mat->nullsp) { MatNullSpaceDestroy(mat->nullsp); }
5746: mat->nullsp = nullsp;
5747: return(0);
5748: }
5752: /*@
5753: MatICCFactor - Performs in-place incomplete Cholesky factorization of matrix.
5755: Collective on Mat
5757: Input Parameters:
5758: + mat - the matrix
5759: . row - row/column permutation
5760: . fill - expected fill factor >= 1.0
5761: - level - level of fill, for ICC(k)
5763: Notes:
5764: Probably really in-place only when level of fill is zero, otherwise allocates
5765: new space to store factored matrix and deletes previous memory.
5767: Most users should employ the simplified KSP interface for linear solvers
5768: instead of working directly with matrix algebra routines such as this.
5769: See, e.g., KSPCreate().
5771: Level: developer
5773: Concepts: matrices^incomplete Cholesky factorization
5774: Concepts: Cholesky factorization
5776: .seealso: MatICCFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
5777: @*/
5778: PetscErrorCode MatICCFactor(Mat mat,IS row,MatFactorInfo* info)
5779: {
5787: if (mat->rmap.N != mat->cmap.N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
5788: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
5789: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
5790: if (!mat->ops->iccfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
5791: MatPreallocated(mat);
5792: (*mat->ops->iccfactor)(mat,row,info);
5793: PetscObjectStateIncrease((PetscObject)mat);
5794: return(0);
5795: }
5799: /*@
5800: MatSetValuesAdic - Sets values computed with ADIC automatic differentiation into a matrix.
5802: Not Collective
5804: Input Parameters:
5805: + mat - the matrix
5806: - v - the values compute with ADIC
5808: Level: developer
5810: Notes:
5811: Must call MatSetColoring() before using this routine. Also this matrix must already
5812: have its nonzero pattern determined.
5814: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
5815: MatSetValues(), MatSetColoring(), MatSetValuesAdifor()
5816: @*/
5817: PetscErrorCode MatSetValuesAdic(Mat mat,void *v)
5818: {
5826: if (!mat->assembled) {
5827: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be already assembled");
5828: }
5830: if (!mat->ops->setvaluesadic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
5831: (*mat->ops->setvaluesadic)(mat,v);
5833: MatView_Private(mat);
5834: PetscObjectStateIncrease((PetscObject)mat);
5835: return(0);
5836: }
5841: /*@
5842: MatSetColoring - Sets a coloring used by calls to MatSetValuesAdic()
5844: Not Collective
5846: Input Parameters:
5847: + mat - the matrix
5848: - coloring - the coloring
5850: Level: developer
5852: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
5853: MatSetValues(), MatSetValuesAdic()
5854: @*/
5855: PetscErrorCode MatSetColoring(Mat mat,ISColoring coloring)
5856: {
5864: if (!mat->assembled) {
5865: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be already assembled");
5866: }
5867: if (!mat->ops->setcoloring) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
5868: (*mat->ops->setcoloring)(mat,coloring);
5869: return(0);
5870: }
5874: /*@
5875: MatSetValuesAdifor - Sets values computed with automatic differentiation into a matrix.
5877: Not Collective
5879: Input Parameters:
5880: + mat - the matrix
5881: . nl - leading dimension of v
5882: - v - the values compute with ADIFOR
5884: Level: developer
5886: Notes:
5887: Must call MatSetColoring() before using this routine. Also this matrix must already
5888: have its nonzero pattern determined.
5890: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
5891: MatSetValues(), MatSetColoring()
5892: @*/
5893: PetscErrorCode MatSetValuesAdifor(Mat mat,PetscInt nl,void *v)
5894: {
5902: if (!mat->assembled) {
5903: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be already assembled");
5904: }
5906: if (!mat->ops->setvaluesadifor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
5907: (*mat->ops->setvaluesadifor)(mat,nl,v);
5909: PetscObjectStateIncrease((PetscObject)mat);
5910: return(0);
5911: }
5915: /*@
5916: MatDiagonalScaleLocal - Scales columns of a matrix given the scaling values including the
5917: ghosted ones.
5919: Not Collective
5921: Input Parameters:
5922: + mat - the matrix
5923: - diag = the diagonal values, including ghost ones
5925: Level: developer
5927: Notes: Works only for MPIAIJ and MPIBAIJ matrices
5928:
5929: .seealso: MatDiagonalScale()
5930: @*/
5931: PetscErrorCode MatDiagonalScaleLocal(Mat mat,Vec diag)
5932: {
5934: PetscMPIInt size;
5941: if (!mat->assembled) {
5942: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be already assembled");
5943: }
5945: MPI_Comm_size(mat->comm,&size);
5946: if (size == 1) {
5947: PetscInt n,m;
5948: VecGetSize(diag,&n);
5949: MatGetSize(mat,0,&m);
5950: if (m == n) {
5951: MatDiagonalScale(mat,0,diag);
5952: } else {
5953: SETERRQ(PETSC_ERR_SUP,"Only supported for sequential matrices when no ghost points/periodic conditions");
5954: }
5955: } else {
5956: PetscErrorCode (*f)(Mat,Vec);
5957: PetscObjectQueryFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",(void (**)(void))&f);
5958: if (f) {
5959: (*f)(mat,diag);
5960: } else {
5961: SETERRQ(PETSC_ERR_SUP,"Only supported for MPIAIJ and MPIBAIJ parallel matrices");
5962: }
5963: }
5965: PetscObjectStateIncrease((PetscObject)mat);
5966: return(0);
5967: }
5971: /*@
5972: MatGetInertia - Gets the inertia from a factored matrix
5974: Collective on Mat
5976: Input Parameter:
5977: . mat - the matrix
5979: Output Parameters:
5980: + nneg - number of negative eigenvalues
5981: . nzero - number of zero eigenvalues
5982: - npos - number of positive eigenvalues
5984: Level: advanced
5986: Notes: Matrix must have been factored by MatCholeskyFactor()
5989: @*/
5990: PetscErrorCode MatGetInertia(Mat mat,PetscInt *nneg,PetscInt *nzero,PetscInt *npos)
5991: {
5997: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
5998: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Numeric factor mat is not assembled");
5999: if (!mat->ops->getinertia) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
6000: (*mat->ops->getinertia)(mat,nneg,nzero,npos);
6001: return(0);
6002: }
6004: /* ----------------------------------------------------------------*/
6007: /*@
6008: MatSolves - Solves A x = b, given a factored matrix, for a collection of vectors
6010: Collective on Mat and Vecs
6012: Input Parameters:
6013: + mat - the factored matrix
6014: - b - the right-hand-side vectors
6016: Output Parameter:
6017: . x - the result vectors
6019: Notes:
6020: The vectors b and x cannot be the same. I.e., one cannot
6021: call MatSolves(A,x,x).
6023: Notes:
6024: Most users should employ the simplified KSP interface for linear solvers
6025: instead of working directly with matrix algebra routines such as this.
6026: See, e.g., KSPCreate().
6028: Level: developer
6030: Concepts: matrices^triangular solves
6032: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd(), MatSolve()
6033: @*/
6034: PetscErrorCode MatSolves(Mat mat,Vecs b,Vecs x)
6035: {
6041: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
6042: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
6043: if (!mat->rmap.N && !mat->cmap.N) return(0);
6045: if (!mat->ops->solves) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
6046: MatPreallocated(mat);
6048: (*mat->ops->solves)(mat,b,x);
6050: return(0);
6051: }
6055: /*@
6056: MatIsSymmetric - Test whether a matrix is symmetric
6058: Collective on Mat
6060: Input Parameter:
6061: + A - the matrix to test
6062: - tol - difference between value and its transpose less than this amount counts as equal (use 0.0 for exact transpose)
6064: Output Parameters:
6065: . flg - the result
6067: Level: intermediate
6069: Concepts: matrix^symmetry
6071: .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsStructurallySymmetric(), MatSetOption(), MatIsSymmetricKnown()
6072: @*/
6073: PetscErrorCode MatIsSymmetric(Mat A,PetscReal tol,PetscTruth *flg)
6074: {
6080: if (!A->symmetric_set) {
6081: if (!A->ops->issymmetric) {
6082: MatType mattype;
6083: MatGetType(A,&mattype);
6084: SETERRQ1(PETSC_ERR_SUP,"Matrix of type <%s> does not support checking for symmetric",mattype);
6085: }
6086: (*A->ops->issymmetric)(A,tol,&A->symmetric);
6087: A->symmetric_set = PETSC_TRUE;
6088: if (A->symmetric) {
6089: A->structurally_symmetric_set = PETSC_TRUE;
6090: A->structurally_symmetric = PETSC_TRUE;
6091: }
6092: }
6093: *flg = A->symmetric;
6094: return(0);
6095: }
6099: /*@
6100: MatIsSymmetricKnown - Checks the flag on the matrix to see if it is symmetric.
6102: Collective on Mat
6104: Input Parameter:
6105: . A - the matrix to check
6107: Output Parameters:
6108: + set - if the symmetric flag is set (this tells you if the next flag is valid)
6109: - flg - the result
6111: Level: advanced
6113: Concepts: matrix^symmetry
6115: Note: Does not check the matrix values directly, so this may return unknown (set = PETSC_FALSE). Use MatIsSymmetric()
6116: if you want it explicitly checked
6118: .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsStructurallySymmetric(), MatSetOption(), MatIsSymmetric()
6119: @*/
6120: PetscErrorCode MatIsSymmetricKnown(Mat A,PetscTruth *set,PetscTruth *flg)
6121: {
6126: if (A->symmetric_set) {
6127: *set = PETSC_TRUE;
6128: *flg = A->symmetric;
6129: } else {
6130: *set = PETSC_FALSE;
6131: }
6132: return(0);
6133: }
6137: /*@
6138: MatIsHermitianKnown - Checks the flag on the matrix to see if it is hermitian.
6140: Collective on Mat
6142: Input Parameter:
6143: . A - the matrix to check
6145: Output Parameters:
6146: + set - if the hermitian flag is set (this tells you if the next flag is valid)
6147: - flg - the result
6149: Level: advanced
6151: Concepts: matrix^symmetry
6153: Note: Does not check the matrix values directly, so this may return unknown (set = PETSC_FALSE). Use MatIsHermitian()
6154: if you want it explicitly checked
6156: .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsStructurallySymmetric(), MatSetOption(), MatIsSymmetric()
6157: @*/
6158: PetscErrorCode MatIsHermitianKnown(Mat A,PetscTruth *set,PetscTruth *flg)
6159: {
6164: if (A->hermitian_set) {
6165: *set = PETSC_TRUE;
6166: *flg = A->hermitian;
6167: } else {
6168: *set = PETSC_FALSE;
6169: }
6170: return(0);
6171: }
6175: /*@
6176: MatIsStructurallySymmetric - Test whether a matrix is structurally symmetric
6178: Collective on Mat
6180: Input Parameter:
6181: . A - the matrix to test
6183: Output Parameters:
6184: . flg - the result
6186: Level: intermediate
6188: Concepts: matrix^symmetry
6190: .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsSymmetric(), MatSetOption()
6191: @*/
6192: PetscErrorCode MatIsStructurallySymmetric(Mat A,PetscTruth *flg)
6193: {
6199: if (!A->structurally_symmetric_set) {
6200: if (!A->ops->isstructurallysymmetric) SETERRQ(PETSC_ERR_SUP,"Matrix does not support checking for structural symmetric");
6201: (*A->ops->isstructurallysymmetric)(A,&A->structurally_symmetric);
6202: A->structurally_symmetric_set = PETSC_TRUE;
6203: }
6204: *flg = A->structurally_symmetric;
6205: return(0);
6206: }
6210: /*@
6211: MatIsHermitian - Test whether a matrix is Hermitian, i.e. it is the complex conjugate of its transpose.
6213: Collective on Mat
6215: Input Parameter:
6216: . A - the matrix to test
6218: Output Parameters:
6219: . flg - the result
6221: Level: intermediate
6223: Concepts: matrix^symmetry
6225: .seealso: MatTranspose(), MatIsTranspose(), MatIsSymmetric(), MatIsStructurallySymmetric(), MatSetOption()
6226: @*/
6227: PetscErrorCode MatIsHermitian(Mat A,PetscTruth *flg)
6228: {
6234: if (!A->hermitian_set) {
6235: if (!A->ops->ishermitian) SETERRQ(PETSC_ERR_SUP,"Matrix does not support checking for being Hermitian");
6236: (*A->ops->ishermitian)(A,&A->hermitian);
6237: A->hermitian_set = PETSC_TRUE;
6238: if (A->hermitian) {
6239: A->structurally_symmetric_set = PETSC_TRUE;
6240: A->structurally_symmetric = PETSC_TRUE;
6241: }
6242: }
6243: *flg = A->hermitian;
6244: return(0);
6245: }
6250: /*@
6251: MatStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
6252: to be communicated to other processors during the MatAssemblyBegin/End() process
6254: Not collective
6256: Input Parameter:
6257: . vec - the vector
6259: Output Parameters:
6260: + nstash - the size of the stash
6261: . reallocs - the number of additional mallocs incurred.
6262: . bnstash - the size of the block stash
6263: - breallocs - the number of additional mallocs incurred.in the block stash
6264:
6265: Level: advanced
6267: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), Mat, MatStashSetInitialSize()
6268:
6269: @*/
6270: PetscErrorCode MatStashGetInfo(Mat mat,PetscInt *nstash,PetscInt *reallocs,PetscInt *bnstash,PetscInt *breallocs)
6271: {
6274: MatStashGetInfo_Private(&mat->stash,nstash,reallocs);
6275: MatStashGetInfo_Private(&mat->bstash,bnstash,breallocs);
6276: return(0);
6277: }
6281: /*@
6282: MatGetVecs - Get vector(s) compatible with the matrix, i.e. with the same
6283: parallel layout
6284:
6285: Collective on Mat
6287: Input Parameter:
6288: . mat - the matrix
6290: Output Parameter:
6291: + right - (optional) vector that the matrix can be multiplied against
6292: - left - (optional) vector that the matrix vector product can be stored in
6294: Level: advanced
6296: .seealso: MatCreate()
6297: @*/
6298: PetscErrorCode MatGetVecs(Mat mat,Vec *right,Vec *left)
6299: {
6305: MatPreallocated(mat);
6306: if (mat->ops->getvecs) {
6307: (*mat->ops->getvecs)(mat,right,left);
6308: } else {
6309: PetscMPIInt size;
6310: MPI_Comm_size(mat->comm, &size);
6311: if (right) {
6312: VecCreate(mat->comm,right);
6313: VecSetSizes(*right,mat->cmap.n,PETSC_DETERMINE);
6314: if (size > 1) {VecSetType(*right,VECMPI);}
6315: else {VecSetType(*right,VECSEQ);}
6316: }
6317: if (left) {
6318: VecCreate(mat->comm,left);
6319: VecSetSizes(*left,mat->rmap.n,PETSC_DETERMINE);
6320: if (size > 1) {VecSetType(*left,VECMPI);}
6321: else {VecSetType(*left,VECSEQ);}
6322: }
6323: }
6324: if (right) {VecSetBlockSize(*right,mat->rmap.bs);}
6325: if (left) {VecSetBlockSize(*left,mat->rmap.bs);}
6326: return(0);
6327: }
6331: /*@
6332: MatFactorInfoInitialize - Initializes a MatFactorInfo data structure
6333: with default values.
6335: Not Collective
6337: Input Parameters:
6338: . info - the MatFactorInfo data structure
6341: Notes: The solvers are generally used through the KSP and PC objects, for example
6342: PCLU, PCILU, PCCHOLESKY, PCICC
6344: Level: developer
6346: .seealso: MatFactorInfo
6347: @*/
6349: PetscErrorCode MatFactorInfoInitialize(MatFactorInfo *info)
6350: {
6354: PetscMemzero(info,sizeof(MatFactorInfo));
6355: return(0);
6356: }
6360: /*@
6361: MatPtAP - Creates the matrix projection C = P^T * A * P
6363: Collective on Mat
6365: Input Parameters:
6366: + A - the matrix
6367: . P - the projection matrix
6368: . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
6369: - fill - expected fill as ratio of nnz(C)/nnz(A)
6371: Output Parameters:
6372: . C - the product matrix
6374: Notes:
6375: C will be created and must be destroyed by the user with MatDestroy().
6377: This routine is currently only implemented for pairs of AIJ matrices and classes
6378: which inherit from AIJ.
6380: Level: intermediate
6382: .seealso: MatPtAPSymbolic(), MatPtAPNumeric(), MatMatMult()
6383: @*/
6384: PetscErrorCode MatPtAP(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
6385: {
6391: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6392: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6395: MatPreallocated(P);
6396: if (!P->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6397: if (P->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6399: if (P->rmap.N!=A->cmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",P->rmap.N,A->cmap.N);
6400: if (fill < 1.0) SETERRQ1(PETSC_ERR_ARG_SIZ,"Expected fill=%G must be >= 1.0",fill);
6401: MatPreallocated(A);
6404: (*A->ops->ptap)(A,P,scall,fill,C);
6407: return(0);
6408: }
6412: /*@
6413: MatPtAPNumeric - Computes the matrix projection C = P^T * A * P
6415: Collective on Mat
6417: Input Parameters:
6418: + A - the matrix
6419: - P - the projection matrix
6421: Output Parameters:
6422: . C - the product matrix
6424: Notes:
6425: C must have been created by calling MatPtAPSymbolic and must be destroyed by
6426: the user using MatDeatroy().
6428: This routine is currently only implemented for pairs of AIJ matrices and classes
6429: which inherit from AIJ. C will be of type MATAIJ.
6431: Level: intermediate
6433: .seealso: MatPtAP(), MatPtAPSymbolic(), MatMatMultNumeric()
6434: @*/
6435: PetscErrorCode MatPtAPNumeric(Mat A,Mat P,Mat C)
6436: {
6442: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6443: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6446: MatPreallocated(P);
6447: if (!P->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6448: if (P->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6451: MatPreallocated(C);
6452: if (C->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6453: if (P->cmap.N!=C->rmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",P->cmap.N,C->rmap.N);
6454: if (P->rmap.N!=A->cmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",P->rmap.N,A->cmap.N);
6455: if (A->rmap.N!=A->cmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix 'A' must be square, %D != %D",A->rmap.N,A->cmap.N);
6456: if (P->cmap.N!=C->cmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",P->cmap.N,C->cmap.N);
6457: MatPreallocated(A);
6460: (*A->ops->ptapnumeric)(A,P,C);
6462: return(0);
6463: }
6467: /*@
6468: MatPtAPSymbolic - Creates the (i,j) structure of the matrix projection C = P^T * A * P
6470: Collective on Mat
6472: Input Parameters:
6473: + A - the matrix
6474: - P - the projection matrix
6476: Output Parameters:
6477: . C - the (i,j) structure of the product matrix
6479: Notes:
6480: C will be created and must be destroyed by the user with MatDestroy().
6482: This routine is currently only implemented for pairs of SeqAIJ matrices and classes
6483: which inherit from SeqAIJ. C will be of type MATSEQAIJ. The product is computed using
6484: this (i,j) structure by calling MatPtAPNumeric().
6486: Level: intermediate
6488: .seealso: MatPtAP(), MatPtAPNumeric(), MatMatMultSymbolic()
6489: @*/
6490: PetscErrorCode MatPtAPSymbolic(Mat A,Mat P,PetscReal fill,Mat *C)
6491: {
6497: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6498: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6499: if (fill <1.0) SETERRQ1(PETSC_ERR_ARG_SIZ,"Expected fill=%G must be >= 1.0",fill);
6502: MatPreallocated(P);
6503: if (!P->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6504: if (P->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6507: if (P->rmap.N!=A->cmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",P->rmap.N,A->cmap.N);
6508: if (A->rmap.N!=A->cmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix 'A' must be square, %D != %D",A->rmap.N,A->cmap.N);
6509: MatPreallocated(A);
6511: (*A->ops->ptapsymbolic)(A,P,fill,C);
6514: MatSetBlockSize(*C,A->rmap.bs);
6516: return(0);
6517: }
6521: /*@
6522: MatMatMult - Performs Matrix-Matrix Multiplication C=A*B.
6524: Collective on Mat
6526: Input Parameters:
6527: + A - the left matrix
6528: . B - the right matrix
6529: . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
6530: - fill - expected fill as ratio of nnz(C)/(nnz(A) + nnz(B))
6532: Output Parameters:
6533: . C - the product matrix
6535: Notes:
6536: C will be created and must be destroyed by the user with MatDestroy().
6537: Unless scall is MAT_REUSE_MATRIX
6539: If you have many matrices with the same non-zero structure to multiply, you
6540: should either
6541: $ 1) use MAT_REUSE_MATRIX in all calls but the first or
6542: $ 2) call MatMatMultSymbolic() once and then MatMatMultNumeric() for each product needed
6544: Level: intermediate
6546: .seealso: MatMatMultSymbolic(), MatMatMultNumeric(), MatPtAP()
6547: @*/
6548: PetscErrorCode MatMatMult(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
6549: {
6551: PetscErrorCode (*fA)(Mat,Mat,MatReuse,PetscReal,Mat*);
6552: PetscErrorCode (*fB)(Mat,Mat,MatReuse,PetscReal,Mat*);
6553: PetscErrorCode (*mult)(Mat,Mat,MatReuse,PetscReal,Mat *)=PETSC_NULL;
6558: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6559: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6562: MatPreallocated(B);
6563: if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6564: if (B->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6566: if (B->rmap.N!=A->cmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap.N,A->cmap.N);
6567: if (fill == PETSC_DEFAULT) fill = 2.0;
6568: if (fill < 1.0) SETERRQ1(PETSC_ERR_ARG_SIZ,"Expected fill=%G must be >= 1.0",fill);
6569: MatPreallocated(A);
6571: fA = A->ops->matmult;
6572: fB = B->ops->matmult;
6573: if (fB == fA) {
6574: if (!fB) SETERRQ1(PETSC_ERR_SUP,"MatMatMult not supported for B of type %s",B->type_name);
6575: mult = fB;
6576: } else {
6577: /* dispatch based on the type of A and B */
6578: char multname[256];
6579: PetscStrcpy(multname,"MatMatMult_");
6580: PetscStrcat(multname,A->type_name);
6581: PetscStrcat(multname,"_");
6582: PetscStrcat(multname,B->type_name);
6583: PetscStrcat(multname,"_C"); /* e.g., multname = "MatMatMult_aij_dense_C" */
6584: PetscObjectQueryFunction((PetscObject)B,multname,(void (**)(void))&mult);
6585: if (!mult) SETERRQ2(PETSC_ERR_ARG_INCOMP,"MatMatMult requires A, %s, to be compatible with B, %s",A->type_name,B->type_name);
6586: }
6588: (*mult)(A,B,scall,fill,C);
6590: return(0);
6591: }
6595: /*@
6596: MatMatMultSymbolic - Performs construction, preallocation, and computes the ij structure
6597: of the matrix-matrix product C=A*B. Call this routine before calling MatMatMultNumeric().
6599: Collective on Mat
6601: Input Parameters:
6602: + A - the left matrix
6603: . B - the right matrix
6604: - fill - expected fill as ratio of nnz(C)/(nnz(A) + nnz(B))
6606: Output Parameters:
6607: . C - the matrix containing the ij structure of product matrix
6609: Notes:
6610: C will be created and must be destroyed by the user with MatDestroy().
6612: This routine is currently implemented for
6613: - pairs of AIJ matrices and classes which inherit from AIJ, C will be of type MATAIJ.
6614: - pairs of AIJ (A) and Dense (B) matrix, C will be of type MATDENSE.
6616: Level: intermediate
6618: .seealso: MatMatMult(), MatMatMultNumeric()
6619: @*/
6620: PetscErrorCode MatMatMultSymbolic(Mat A,Mat B,PetscReal fill,Mat *C)
6621: {
6623: PetscErrorCode (*Asymbolic)(Mat,Mat,PetscReal,Mat *);
6624: PetscErrorCode (*Bsymbolic)(Mat,Mat,PetscReal,Mat *);
6625: PetscErrorCode (*symbolic)(Mat,Mat,PetscReal,Mat *)=PETSC_NULL;
6630: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6631: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6635: MatPreallocated(B);
6636: if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6637: if (B->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6640: if (B->rmap.N!=A->cmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap.N,A->cmap.N);
6641: if (fill < 1.0) SETERRQ1(PETSC_ERR_ARG_SIZ,"Expected fill=%G must be > 1.0",fill);
6642: MatPreallocated(A);
6643:
6644: Asymbolic = A->ops->matmultsymbolic;
6645: Bsymbolic = B->ops->matmultsymbolic;
6646: if (Asymbolic == Bsymbolic){
6647: if (!Bsymbolic) SETERRQ1(PETSC_ERR_SUP,"C=A*B not implemented for B of type %s",B->type_name);
6648: symbolic = Bsymbolic;
6649: } else { /* dispatch based on the type of A and B */
6650: char symbolicname[256];
6651: PetscStrcpy(symbolicname,"MatMatMultSymbolic_");
6652: PetscStrcat(symbolicname,A->type_name);
6653: PetscStrcat(symbolicname,"_");
6654: PetscStrcat(symbolicname,B->type_name);
6655: PetscStrcat(symbolicname,"_C");
6656: PetscObjectQueryFunction((PetscObject)B,symbolicname,(void (**)(void))&symbolic);
6657: if (!symbolic) SETERRQ2(PETSC_ERR_ARG_INCOMP,"MatMatMultSymbolic requires A, %s, to be compatible with B, %s",A->type_name,B->type_name);
6658: }
6660: (*symbolic)(A,B,fill,C);
6662: return(0);
6663: }
6667: /*@
6668: MatMatMultNumeric - Performs the numeric matrix-matrix product.
6669: Call this routine after first calling MatMatMultSymbolic().
6671: Collective on Mat
6673: Input Parameters:
6674: + A - the left matrix
6675: - B - the right matrix
6677: Output Parameters:
6678: . C - the product matrix, whose ij structure was defined from MatMatMultSymbolic().
6680: Notes:
6681: C must have been created with MatMatMultSymbolic.
6683: This routine is currently implemented for
6684: - pairs of AIJ matrices and classes which inherit from AIJ, C will be of type MATAIJ.
6685: - pairs of AIJ (A) and Dense (B) matrix, C will be of type MATDENSE.
6687: Level: intermediate
6689: .seealso: MatMatMult(), MatMatMultSymbolic()
6690: @*/
6691: PetscErrorCode MatMatMultNumeric(Mat A,Mat B,Mat C)
6692: {
6694: PetscErrorCode (*Anumeric)(Mat,Mat,Mat);
6695: PetscErrorCode (*Bnumeric)(Mat,Mat,Mat);
6696: PetscErrorCode (*numeric)(Mat,Mat,Mat)=PETSC_NULL;
6701: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6702: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6706: MatPreallocated(B);
6707: if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6708: if (B->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6712: MatPreallocated(C);
6713: if (!C->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6714: if (C->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6716: if (B->cmap.N!=C->cmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->cmap.N,C->cmap.N);
6717: if (B->rmap.N!=A->cmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap.N,A->cmap.N);
6718: if (A->rmap.N!=C->rmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",A->rmap.N,C->rmap.N);
6719: MatPreallocated(A);
6721: Anumeric = A->ops->matmultnumeric;
6722: Bnumeric = B->ops->matmultnumeric;
6723: if (Anumeric == Bnumeric){
6724: if (!Bnumeric) SETERRQ1(PETSC_ERR_SUP,"MatMatMultNumeric not supported for B of type %s",B->type_name);
6725: numeric = Bnumeric;
6726: } else {
6727: char numericname[256];
6728: PetscStrcpy(numericname,"MatMatMultNumeric_");
6729: PetscStrcat(numericname,A->type_name);
6730: PetscStrcat(numericname,"_");
6731: PetscStrcat(numericname,B->type_name);
6732: PetscStrcat(numericname,"_C");
6733: PetscObjectQueryFunction((PetscObject)B,numericname,(void (**)(void))&numeric);
6734: if (!numeric)
6735: SETERRQ2(PETSC_ERR_ARG_INCOMP,"MatMatMultNumeric requires A, %s, to be compatible with B, %s",A->type_name,B->type_name);
6736: }
6738: (*numeric)(A,B,C);
6740: return(0);
6741: }
6745: /*@
6746: MatMatMultTranspose - Performs Matrix-Matrix Multiplication C=A^T*B.
6748: Collective on Mat
6750: Input Parameters:
6751: + A - the left matrix
6752: . B - the right matrix
6753: . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
6754: - fill - expected fill as ratio of nnz(C)/(nnz(A) + nnz(B))
6756: Output Parameters:
6757: . C - the product matrix
6759: Notes:
6760: C will be created and must be destroyed by the user with MatDestroy().
6762: This routine is currently only implemented for pairs of SeqAIJ matrices and classes
6763: which inherit from SeqAIJ. C will be of type MATSEQAIJ.
6765: Level: intermediate
6767: .seealso: MatMatMultTransposeSymbolic(), MatMatMultTransposeNumeric(), MatPtAP()
6768: @*/
6769: PetscErrorCode MatMatMultTranspose(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
6770: {
6772: PetscErrorCode (*fA)(Mat,Mat,MatReuse,PetscReal,Mat*);
6773: PetscErrorCode (*fB)(Mat,Mat,MatReuse,PetscReal,Mat*);
6778: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6779: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6782: MatPreallocated(B);
6783: if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
6784: if (B->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
6786: if (B->rmap.N!=A->rmap.N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap.N,A->rmap.N);
6787: if (fill < 1.0) SETERRQ1(PETSC_ERR_ARG_SIZ,"Expected fill=%G must be > 1.0",fill);
6788: MatPreallocated(A);
6790: fA = A->ops->matmulttranspose;
6791: if (!fA) SETERRQ1(PETSC_ERR_SUP,"MatMatMultTranspose not supported for A of type %s",A->type_name);
6792: fB = B->ops->matmulttranspose;
6793: if (!fB) SETERRQ1(PETSC_ERR_SUP,"MatMatMultTranspose not supported for B of type %s",B->type_name);
6794: if (fB!=fA) SETERRQ2(PETSC_ERR_ARG_INCOMP,"MatMatMultTranspose requires A, %s, to be compatible with B, %s",A->type_name,B->type_name);
6797: (*A->ops->matmulttranspose)(A,B,scall,fill,C);
6799:
6800: return(0);
6801: }