Actual source code: matrix.c
1: /*$Id: matrix.c,v 1.414 2001/09/28 18:57:28 balay Exp $*/
3: /*
4: This is where the abstract matrix operations are defined
5: */
7: #include src/mat/matimpl.h
8: #include src/vec/vecimpl.h
10: /* Logging support */
11: int MAT_COOKIE;
12: int MATSNESMFCTX_COOKIE;
13: int MAT_Mult, MAT_MultMatrixFree, MAT_Mults, MAT_MultConstrained, MAT_MultAdd, MAT_MultTranspose;
14: int MAT_MultTransposeConstrained, MAT_MultTransposeAdd, MAT_Solve, MAT_Solves, MAT_SolveAdd, MAT_SolveTranspose;
15: int MAT_SolveTransposeAdd, MAT_Relax, MAT_ForwardSolve, MAT_BackwardSolve, MAT_LUFactor, MAT_LUFactorSymbolic;
16: int MAT_LUFactorNumeric, MAT_CholeskyFactor, MAT_CholeskyFactorSymbolic, MAT_CholeskyFactorNumeric, MAT_ILUFactor;
17: int MAT_ILUFactorSymbolic, MAT_ICCFactorSymbolic, MAT_Copy, MAT_Convert, MAT_Scale, MAT_AssemblyBegin;
18: int MAT_AssemblyEnd, MAT_SetValues, MAT_GetValues, MAT_GetRow, MAT_GetSubMatrices, MAT_GetColoring, MAT_GetOrdering;
19: int MAT_IncreaseOverlap, MAT_Partitioning, MAT_ZeroEntries, MAT_Load, MAT_View, MAT_AXPY, MAT_FDColoringCreate;
20: int MAT_FDColoringApply,MAT_Transpose;
22: #undef __FUNCT__
24: /*@C
25: MatGetRow - Gets a row of a matrix. You MUST call MatRestoreRow()
26: for each row that you get to ensure that your application does
27: not bleed memory.
29: Not Collective
31: Input Parameters:
32: + mat - the matrix
33: - row - the row to get
35: Output Parameters:
36: + ncols - the number of nonzeros in the row
37: . cols - if not NULL, the column numbers
38: - vals - if not NULL, the values
40: Notes:
41: This routine is provided for people who need to have direct access
42: to the structure of a matrix. We hope that we provide enough
43: high-level matrix routines that few users will need it.
45: MatGetRow() always returns 0-based column indices, regardless of
46: whether the internal representation is 0-based (default) or 1-based.
48: For better efficiency, set cols and/or vals to PETSC_NULL if you do
49: not wish to extract these quantities.
51: The user can only examine the values extracted with MatGetRow();
52: the values cannot be altered. To change the matrix entries, one
53: must use MatSetValues().
55: You can only have one call to MatGetRow() outstanding for a particular
56: matrix at a time, per processor. MatGetRow() can only obtained rows
57: associated with the given processor, it cannot get rows from the
58: other processors; for that we suggest using MatGetSubMatrices(), then
59: MatGetRow() on the submatrix. The row indix passed to MatGetRows()
60: is in the global number of rows.
62: Fortran Notes:
63: The calling sequence from Fortran is
64: .vb
65: MatGetRow(matrix,row,ncols,cols,values,ierr)
66: Mat matrix (input)
67: integer row (input)
68: integer ncols (output)
69: integer cols(maxcols) (output)
70: double precision (or double complex) values(maxcols) output
71: .ve
72: where maxcols >= maximum nonzeros in any row of the matrix.
74: Caution:
75: Do not try to change the contents of the output arrays (cols and vals).
76: In some cases, this may corrupt the matrix.
78: Level: advanced
80: Concepts: matrices^row access
82: .seealso: MatRestoreRow(), MatSetValues(), MatGetValues(), MatGetSubmatrices(), MatGetDiagonal()
83: @*/
84: int MatGetRow(Mat mat,int row,int *ncols,int **cols,PetscScalar **vals)
85: {
86: int ierr;
92: MatPreallocated(mat);
93: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
94: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
95: if (!mat->ops->getrow) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
96: PetscLogEventBegin(MAT_GetRow,mat,0,0,0);
97: (*mat->ops->getrow)(mat,row,ncols,cols,vals);
98: PetscLogEventEnd(MAT_GetRow,mat,0,0,0);
99: return(0);
100: }
102: #undef __FUNCT__
104: /*@C
105: MatRestoreRow - Frees any temporary space allocated by MatGetRow().
107: Not Collective
109: Input Parameters:
110: + mat - the matrix
111: . row - the row to get
112: . ncols, cols - the number of nonzeros and their columns
113: - vals - if nonzero the column values
115: Notes:
116: This routine should be called after you have finished examining the entries.
118: Fortran Notes:
119: The calling sequence from Fortran is
120: .vb
121: MatRestoreRow(matrix,row,ncols,cols,values,ierr)
122: Mat matrix (input)
123: integer row (input)
124: integer ncols (output)
125: integer cols(maxcols) (output)
126: double precision (or double complex) values(maxcols) output
127: .ve
128: Where maxcols >= maximum nonzeros in any row of the matrix.
130: In Fortran MatRestoreRow() MUST be called after MatGetRow()
131: before another call to MatGetRow() can be made.
133: Level: advanced
135: .seealso: MatGetRow()
136: @*/
137: int MatRestoreRow(Mat mat,int row,int *ncols,int **cols,PetscScalar **vals)
138: {
144: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
145: if (!mat->ops->restorerow) return(0);
146: (*mat->ops->restorerow)(mat,row,ncols,cols,vals);
147: return(0);
148: }
150: #undef __FUNCT__
152: /*@C
153: MatView - Visualizes a matrix object.
155: Collective on Mat
157: Input Parameters:
158: + mat - the matrix
159: - ptr - visualization context
161: Notes:
162: The available visualization contexts include
163: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
164: . PETSC_VIEWER_STDOUT_WORLD - synchronized standard
165: output where only the first processor opens
166: the file. All other processors send their
167: data to the first processor to print.
168: - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure
170: The user can open alternative visualization contexts with
171: + PetscViewerASCIIOpen() - Outputs matrix to a specified file
172: . PetscViewerBinaryOpen() - Outputs matrix in binary to a
173: specified file; corresponding input uses MatLoad()
174: . PetscViewerDrawOpen() - Outputs nonzero matrix structure to
175: an X window display
176: - PetscViewerSocketOpen() - Outputs matrix to Socket viewer.
177: Currently only the sequential dense and AIJ
178: matrix types support the Socket viewer.
180: The user can call PetscViewerSetFormat() to specify the output
181: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
182: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
183: + PETSC_VIEWER_ASCII_DEFAULT - default, prints matrix contents
184: . PETSC_VIEWER_ASCII_MATLAB - prints matrix contents in Matlab format
185: . PETSC_VIEWER_ASCII_DENSE - prints entire matrix including zeros
186: . PETSC_VIEWER_ASCII_COMMON - prints matrix contents, using a sparse
187: format common among all matrix types
188: . PETSC_VIEWER_ASCII_IMPL - prints matrix contents, using an implementation-specific
189: format (which is in many cases the same as the default)
190: . PETSC_VIEWER_ASCII_INFO - prints basic information about the matrix
191: size and structure (not the matrix entries)
192: . PETSC_VIEWER_ASCII_INFO_DETAIL - prints more detailed information about
193: the matrix structure
195: Options Database Keys:
196: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
197: . -mat_view_info_detailed - Prints more detailed info
198: . -mat_view - Prints matrix in ASCII format
199: . -mat_view_matlab - Prints matrix in Matlab format
200: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
201: . -display <name> - Sets display name (default is host)
202: . -draw_pause <sec> - Sets number of seconds to pause after display
203: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
204: . -viewer_socket_machine <machine>
205: . -viewer_socket_port <port>
206: . -mat_view_binary - save matrix to file in binary format
207: - -viewer_binary_filename <name>
208: Level: beginner
210: Concepts: matrices^viewing
211: Concepts: matrices^plotting
212: Concepts: matrices^printing
214: .seealso: PetscViewerSetFormat(), PetscViewerASCIIOpen(), PetscViewerDrawOpen(),
215: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), MatLoad()
216: @*/
217: int MatView(Mat mat,PetscViewer viewer)
218: {
219: int ierr,rows,cols;
220: PetscTruth isascii;
221: char *cstr;
222: PetscViewerFormat format;
227: MatPreallocated(mat);
228: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(mat->comm);
231: if (!mat->assembled) SETERRQ(1,"Must call MatAssemblyBegin/End() before viewing matrix");
233: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
234: if (isascii) {
235: PetscViewerGetFormat(viewer,&format);
236: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
237: if (mat->prefix) {
238: PetscViewerASCIIPrintf(viewer,"Matrix Object:(%s)n",mat->prefix);
239: } else {
240: PetscViewerASCIIPrintf(viewer,"Matrix Object:n");
241: }
242: PetscViewerASCIIPushTab(viewer);
243: MatGetType(mat,&cstr);
244: MatGetSize(mat,&rows,&cols);
245: PetscViewerASCIIPrintf(viewer,"type=%s, rows=%d, cols=%dn",cstr,rows,cols);
246: if (mat->ops->getinfo) {
247: MatInfo info;
248: MatGetInfo(mat,MAT_GLOBAL_SUM,&info);
249: PetscViewerASCIIPrintf(viewer,"total: nonzeros=%d, allocated nonzeros=%dn",
250: (int)info.nz_used,(int)info.nz_allocated);
251: }
252: }
253: }
254: if (mat->ops->view) {
255: PetscViewerASCIIPushTab(viewer);
256: (*mat->ops->view)(mat,viewer);
257: PetscViewerASCIIPopTab(viewer);
258: } else if (!isascii) {
259: SETERRQ1(1,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
260: }
261: if (isascii) {
262: PetscViewerGetFormat(viewer,&format);
263: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
264: PetscViewerASCIIPopTab(viewer);
265: }
266: }
267: return(0);
268: }
270: #undef __FUNCT__
272: /*@C
273: MatScaleSystem - Scale a vector solution and right hand side to
274: match the scaling of a scaled matrix.
275:
276: Collective on Mat
278: Input Parameter:
279: + mat - the matrix
280: . x - solution vector (or PETSC_NULL)
281: - b - right hand side vector (or PETSC_NULL)
284: Notes:
285: For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
286: internally scaled, so this does nothing. For MPIROWBS it
287: permutes and diagonally scales.
289: The SLES methods automatically call this routine when required
290: (via PCPreSolve()) so it is rarely used directly.
292: Level: Developer
294: Concepts: matrices^scaling
296: .seealso: MatUseScaledForm(), MatUnScaleSystem()
297: @*/
298: int MatScaleSystem(Mat mat,Vec x,Vec b)
299: {
305: MatPreallocated(mat);
309: if (mat->ops->scalesystem) {
310: (*mat->ops->scalesystem)(mat,x,b);
311: }
312: return(0);
313: }
315: #undef __FUNCT__
317: /*@C
318: MatUnScaleSystem - Unscales a vector solution and right hand side to
319: match the original scaling of a scaled matrix.
320:
321: Collective on Mat
323: Input Parameter:
324: + mat - the matrix
325: . x - solution vector (or PETSC_NULL)
326: - b - right hand side vector (or PETSC_NULL)
329: Notes:
330: For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
331: internally scaled, so this does nothing. For MPIROWBS it
332: permutes and diagonally scales.
334: The SLES methods automatically call this routine when required
335: (via PCPreSolve()) so it is rarely used directly.
337: Level: Developer
339: .seealso: MatUseScaledForm(), MatScaleSystem()
340: @*/
341: int MatUnScaleSystem(Mat mat,Vec x,Vec b)
342: {
348: MatPreallocated(mat);
351: if (mat->ops->unscalesystem) {
352: (*mat->ops->unscalesystem)(mat,x,b);
353: }
354: return(0);
355: }
357: #undef __FUNCT__
359: /*@C
360: MatUseScaledForm - For matrix storage formats that scale the
361: matrix (for example MPIRowBS matrices are diagonally scaled on
362: assembly) indicates matrix operations (MatMult() etc) are
363: applied using the scaled matrix.
364:
365: Collective on Mat
367: Input Parameter:
368: + mat - the matrix
369: - scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for
370: applying the original matrix
372: Notes:
373: For scaled matrix formats, applying the original, unscaled matrix
374: will be slightly more expensive
376: Level: Developer
378: .seealso: MatScaleSystem(), MatUnScaleSystem()
379: @*/
380: int MatUseScaledForm(Mat mat,PetscTruth scaled)
381: {
387: MatPreallocated(mat);
388: if (mat->ops->usescaledform) {
389: (*mat->ops->usescaledform)(mat,scaled);
390: }
391: return(0);
392: }
394: #undef __FUNCT__
396: /*@C
397: MatDestroy - Frees space taken by a matrix.
398:
399: Collective on Mat
401: Input Parameter:
402: . A - the matrix
404: Level: beginner
406: @*/
407: int MatDestroy(Mat A)
408: {
414: MatPreallocated(A);
415: if (--A->refct > 0) return(0);
417: /* if memory was published with AMS then destroy it */
418: PetscObjectDepublish(A);
419: if (A->mapping) {
420: ISLocalToGlobalMappingDestroy(A->mapping);
421: }
422: if (A->bmapping) {
423: ISLocalToGlobalMappingDestroy(A->bmapping);
424: }
425: if (A->rmap) {
426: PetscMapDestroy(A->rmap);
427: }
428: if (A->cmap) {
429: PetscMapDestroy(A->cmap);
430: }
432: (*A->ops->destroy)(A);
433: PetscLogObjectDestroy(A);
434: PetscHeaderDestroy(A);
435: return(0);
436: }
438: #undef __FUNCT__
440: /*@
441: MatValid - Checks whether a matrix object is valid.
443: Collective on Mat
445: Input Parameter:
446: . m - the matrix to check
448: Output Parameter:
449: flg - flag indicating matrix status, either
450: PETSC_TRUE if matrix is valid, or PETSC_FALSE otherwise.
452: Level: developer
454: Concepts: matrices^validity
455: @*/
456: int MatValid(Mat m,PetscTruth *flg)
457: {
460: if (!m) *flg = PETSC_FALSE;
461: else if (m->cookie != MAT_COOKIE) *flg = PETSC_FALSE;
462: else *flg = PETSC_TRUE;
463: return(0);
464: }
466: #undef __FUNCT__
468: /*@
469: MatSetValues - Inserts or adds a block of values into a matrix.
470: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
471: MUST be called after all calls to MatSetValues() have been completed.
473: Not Collective
475: Input Parameters:
476: + mat - the matrix
477: . v - a logically two-dimensional array of values
478: . m, idxm - the number of rows and their global indices
479: . n, idxn - the number of columns and their global indices
480: - addv - either ADD_VALUES or INSERT_VALUES, where
481: ADD_VALUES adds values to any existing entries, and
482: INSERT_VALUES replaces existing entries with new values
484: Notes:
485: By default the values, v, are row-oriented and unsorted.
486: See MatSetOption() for other options.
488: Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES
489: options cannot be mixed without intervening calls to the assembly
490: routines.
492: MatSetValues() uses 0-based row and column numbers in Fortran
493: as well as in C.
495: Negative indices may be passed in idxm and idxn, these rows and columns are
496: simply ignored. This allows easily inserting element stiffness matrices
497: with homogeneous Dirchlet boundary conditions that you don't want represented
498: in the matrix.
500: Efficiency Alert:
501: The routine MatSetValuesBlocked() may offer much better efficiency
502: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
504: Level: beginner
506: Concepts: matrices^putting entries in
508: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
509: @*/
510: int MatSetValues(Mat mat,int m,int *idxm,int n,int *idxn,PetscScalar *v,InsertMode addv)
511: {
515: if (!m || !n) return(0); /* no values to insert */
518: MatPreallocated(mat);
522: if (mat->insertmode == NOT_SET_VALUES) {
523: mat->insertmode = addv;
524: }
525: #if defined(PETSC_USE_BOPT_g)
526: else if (mat->insertmode != addv) {
527: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
528: }
529: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
530: #endif
532: if (mat->assembled) {
533: mat->was_assembled = PETSC_TRUE;
534: mat->assembled = PETSC_FALSE;
535: }
536: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
537: if (!mat->ops->setvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
538: (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);
539: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
540: return(0);
541: }
543: #undef __FUNCT__
545: /*@C
546: MatSetValuesStencil - Inserts or adds a block of values into a matrix.
547: Using structured grid indexing
549: Not Collective
551: Input Parameters:
552: + mat - the matrix
553: . v - a logically two-dimensional array of values
554: . m - number of rows being entered
555: . idxm - grid coordinates (and component number when dof > 1) for matrix rows being entered
556: . n - number of columns being entered
557: . idxn - grid coordinates (and component number when dof > 1) for matrix columns being entered
558: - addv - either ADD_VALUES or INSERT_VALUES, where
559: ADD_VALUES adds values to any existing entries, and
560: INSERT_VALUES replaces existing entries with new values
562: Notes:
563: By default the values, v, are row-oriented and unsorted.
564: See MatSetOption() for other options.
566: Calls to MatSetValuesStencil() with the INSERT_VALUES and ADD_VALUES
567: options cannot be mixed without intervening calls to the assembly
568: routines.
570: The grid coordinates are across the entire grid, not just the local portion
572: MatSetValuesStencil() uses 0-based row and column numbers in Fortran
573: as well as in C.
575: In order to use this routine you must either obtain the matrix with DAGetMatrix()
576: or call MatSetLocalToGlobalMapping() and MatSetStencil() first.
578: In Fortran idxm and idxn should be declared as
579: $ MatStencil idxm(4,m),idxn(4,n)
580: and the values inserted using
581: $ idxm(MatStencil_i,1) = i
582: $ idxm(MatStencil_j,1) = j
583: $ idxm(MatStencil_k,1) = k
584: $ idxm(MatStencil_c,1) = c
585: etc
587: Negative indices may be passed in idxm and idxn, these rows and columns are
588: simply ignored. This allows easily inserting element stiffness matrices
589: with homogeneous Dirchlet boundary conditions that you don't want represented
590: in the matrix.
592: Inspired by the structured grid interface to the HYPRE package
593: (http://www.llnl.gov/CASC/hypre)
595: Efficiency Alert:
596: The routine MatSetValuesBlockedStencil() may offer much better efficiency
597: for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
599: Level: beginner
601: Concepts: matrices^putting entries in
603: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
604: MatSetValues(), MatSetValuesBlockedStencil(), MatSetStencil(), DAGetMatrix()
605: @*/
606: int MatSetValuesStencil(Mat mat,int m,MatStencil *idxm,int n,MatStencil *idxn,PetscScalar *v,InsertMode addv)
607: {
608: int j,i,ierr,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
609: int *starts = mat->stencil.starts,*dxm = (int*)idxm,*dxn = (int*)idxn,sdim = dim - (1 - (int)mat->stencil.noc);
612: if (!m || !n) return(0); /* no values to insert */
619: if (m > 128) SETERRQ1(1,"Can only set 128 rows at a time; trying to set %d",m);
620: if (n > 128) SETERRQ1(1,"Can only set 256 columns at a time; trying to set %d",n);
622: for (i=0; i<m; i++) {
623: for (j=0; j<3-sdim; j++) dxm++;
624: tmp = *dxm++ - starts[0];
625: for (j=0; j<dim-1; j++) {
626: tmp = tmp*dims[j] + *dxm++ - starts[j+1];
627: }
628: if (mat->stencil.noc) dxm++;
629: jdxm[i] = tmp;
630: }
631: for (i=0; i<n; i++) {
632: for (j=0; j<3-sdim; j++) dxn++;
633: tmp = *dxn++ - starts[0];
634: for (j=0; j<dim-1; j++) {
635: tmp = tmp*dims[j] + *dxn++ - starts[j+1];
636: }
637: if (mat->stencil.noc) dxn++;
638: jdxn[i] = tmp;
639: }
640: MatSetValuesLocal(mat,m,jdxm,n,jdxn,v,addv);
641: return(0);
642: }
644: #undef __FUNCT__
646: /*@
647: MatSetStencil - Sets the grid information for setting values into a matrix via
648: MatSetStencil()
650: Not Collective
652: Input Parameters:
653: + mat - the matrix
654: . dim - dimension of the grid 1,2, or 3
655: . dims - number of grid points in x, y, and z direction, including ghost points on your processor
656: . starts - starting point of ghost nodes on your processor in x, y, and z direction
657: - dof - number of degrees of freedom per node
660: Inspired by the structured grid interface to the HYPRE package
661: (www.llnl.gov/CASC/hyper)
663: Level: beginner
665: Concepts: matrices^putting entries in
667: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
668: MatSetValues(), MatSetValuesBlockedStencil(), MatSetValuesStencil()
669: @*/
670: int MatSetStencil(Mat mat,int dim,int *dims,int *starts,int dof)
671: {
672: int i;
679: mat->stencil.dim = dim + (dof > 1);
680: for (i=0; i<dim; i++) {
681: mat->stencil.dims[i] = dims[dim-i-1]; /* copy the values in backwards */
682: mat->stencil.starts[i] = starts[dim-i-1];
683: }
684: mat->stencil.dims[dim] = dof;
685: mat->stencil.starts[dim] = 0;
686: mat->stencil.noc = (PetscTruth)(dof == 1);
687: return(0);
688: }
690: #undef __FUNCT__
692: /*@
693: MatSetValuesBlocked - Inserts or adds a block of values into a matrix.
695: Not Collective
697: Input Parameters:
698: + mat - the matrix
699: . v - a logically two-dimensional array of values
700: . m, idxm - the number of block rows and their global block indices
701: . n, idxn - the number of block columns and their global block indices
702: - addv - either ADD_VALUES or INSERT_VALUES, where
703: ADD_VALUES adds values to any existing entries, and
704: INSERT_VALUES replaces existing entries with new values
706: Notes:
707: By default the values, v, are row-oriented and unsorted. So the layout of
708: v is the same as for MatSetValues(). See MatSetOption() for other options.
710: Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
711: options cannot be mixed without intervening calls to the assembly
712: routines.
714: MatSetValuesBlocked() uses 0-based row and column numbers in Fortran
715: as well as in C.
717: Negative indices may be passed in idxm and idxn, these rows and columns are
718: simply ignored. This allows easily inserting element stiffness matrices
719: with homogeneous Dirchlet boundary conditions that you don't want represented
720: in the matrix.
722: Each time an entry is set within a sparse matrix via MatSetValues(),
723: internal searching must be done to determine where to place the the
724: data in the matrix storage space. By instead inserting blocks of
725: entries via MatSetValuesBlocked(), the overhead of matrix assembly is
726: reduced.
728: Restrictions:
729: MatSetValuesBlocked() is currently supported only for the block AIJ
730: matrix format (MATSEQBAIJ and MATMPIBAIJ, which are created via
731: MatCreateSeqBAIJ() and MatCreateMPIBAIJ()).
733: Level: intermediate
735: Concepts: matrices^putting entries in blocked
737: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
738: @*/
739: int MatSetValuesBlocked(Mat mat,int m,int *idxm,int n,int *idxn,PetscScalar *v,InsertMode addv)
740: {
744: if (!m || !n) return(0); /* no values to insert */
747: MatPreallocated(mat);
751: if (mat->insertmode == NOT_SET_VALUES) {
752: mat->insertmode = addv;
753: }
754: #if defined(PETSC_USE_BOPT_g)
755: else if (mat->insertmode != addv) {
756: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
757: }
758: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
759: #endif
761: if (mat->assembled) {
762: mat->was_assembled = PETSC_TRUE;
763: mat->assembled = PETSC_FALSE;
764: }
765: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
766: if (!mat->ops->setvaluesblocked) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
767: (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);
768: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
769: return(0);
770: }
772: /*MC
773: MatSetValue - Set a single entry into a matrix.
775: Synopsis:
776: int MatSetValue(Mat m,int row,int col,PetscScalar value,InsertMode mode);
778: Not collective
780: Input Parameters:
781: + m - the matrix
782: . row - the row location of the entry
783: . col - the column location of the entry
784: . value - the value to insert
785: - mode - either INSERT_VALUES or ADD_VALUES
787: Notes:
788: For efficiency one should use MatSetValues() and set several or many
789: values simultaneously if possible.
791: Note that VecSetValue() does NOT return an error code (since this
792: is checked internally).
794: Level: beginner
796: .seealso: MatSetValues()
797: M*/
799: #undef __FUNCT__
801: /*@
802: MatGetValues - Gets a block of values from a matrix.
804: Not Collective; currently only returns a local block
806: Input Parameters:
807: + mat - the matrix
808: . v - a logically two-dimensional array for storing the values
809: . m, idxm - the number of rows and their global indices
810: - n, idxn - the number of columns and their global indices
812: Notes:
813: The user must allocate space (m*n PetscScalars) for the values, v.
814: The values, v, are then returned in a row-oriented format,
815: analogous to that used by default in MatSetValues().
817: MatGetValues() uses 0-based row and column numbers in
818: Fortran as well as in C.
820: MatGetValues() requires that the matrix has been assembled
821: with MatAssemblyBegin()/MatAssemblyEnd(). Thus, calls to
822: MatSetValues() and MatGetValues() CANNOT be made in succession
823: without intermediate matrix assembly.
825: Level: advanced
827: Concepts: matrices^accessing values
829: .seealso: MatGetRow(), MatGetSubmatrices(), MatSetValues()
830: @*/
831: int MatGetValues(Mat mat,int m,int *idxm,int n,int *idxn,PetscScalar *v)
832: {
838: MatPreallocated(mat);
842: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
843: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
844: if (!mat->ops->getvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
846: PetscLogEventBegin(MAT_GetValues,mat,0,0,0);
847: (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);
848: PetscLogEventEnd(MAT_GetValues,mat,0,0,0);
849: return(0);
850: }
852: #undef __FUNCT__
854: /*@
855: MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
856: the routine MatSetValuesLocal() to allow users to insert matrix entries
857: using a local (per-processor) numbering.
859: Not Collective
861: Input Parameters:
862: + x - the matrix
863: - mapping - mapping created with ISLocalToGlobalMappingCreate()
864: or ISLocalToGlobalMappingCreateIS()
866: Level: intermediate
868: Concepts: matrices^local to global mapping
869: Concepts: local to global mapping^for matrices
871: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
872: @*/
873: int MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping)
874: {
879: MatPreallocated(x);
881: if (x->mapping) {
882: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
883: }
885: if (x->ops->setlocaltoglobalmapping) {
886: (*x->ops->setlocaltoglobalmapping)(x,mapping);
887: } else {
888: x->mapping = mapping;
889: PetscObjectReference((PetscObject)mapping);
890: }
891: return(0);
892: }
894: #undef __FUNCT__
896: /*@
897: MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use
898: by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
899: entries using a local (per-processor) numbering.
901: Not Collective
903: Input Parameters:
904: + x - the matrix
905: - mapping - mapping created with ISLocalToGlobalMappingCreate() or
906: ISLocalToGlobalMappingCreateIS()
908: Level: intermediate
910: Concepts: matrices^local to global mapping blocked
911: Concepts: local to global mapping^for matrices, blocked
913: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
914: MatSetValuesBlocked(), MatSetValuesLocal()
915: @*/
916: int MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping mapping)
917: {
922: MatPreallocated(x);
924: if (x->bmapping) {
925: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
926: }
927:
928: x->bmapping = mapping;
929: PetscObjectReference((PetscObject)mapping);
930: return(0);
931: }
933: #undef __FUNCT__
935: /*@
936: MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
937: using a local ordering of the nodes.
939: Not Collective
941: Input Parameters:
942: + x - the matrix
943: . nrow, irow - number of rows and their local indices
944: . ncol, icol - number of columns and their local indices
945: . y - a logically two-dimensional array of values
946: - addv - either INSERT_VALUES or ADD_VALUES, where
947: ADD_VALUES adds values to any existing entries, and
948: INSERT_VALUES replaces existing entries with new values
950: Notes:
951: Before calling MatSetValuesLocal(), the user must first set the
952: local-to-global mapping by calling MatSetLocalToGlobalMapping().
954: Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES
955: options cannot be mixed without intervening calls to the assembly
956: routines.
958: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
959: MUST be called after all calls to MatSetValuesLocal() have been completed.
961: Level: intermediate
963: Concepts: matrices^putting entries in with local numbering
965: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping(),
966: MatSetValueLocal()
967: @*/
968: int MatSetValuesLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,PetscScalar *y,InsertMode addv)
969: {
970: int ierr,irowm[2048],icolm[2048];
975: MatPreallocated(mat);
980: if (mat->insertmode == NOT_SET_VALUES) {
981: mat->insertmode = addv;
982: }
983: #if defined(PETSC_USE_BOPT_g)
984: else if (mat->insertmode != addv) {
985: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
986: }
987: if (!mat->ops->setvalueslocal && (nrow > 2048 || ncol > 2048)) {
988: SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
989: }
990: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
991: #endif
993: if (mat->assembled) {
994: mat->was_assembled = PETSC_TRUE;
995: mat->assembled = PETSC_FALSE;
996: }
997: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
998: if (!mat->ops->setvalueslocal) {
999: ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm);
1000: ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);
1001: (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);
1002: } else {
1003: (*mat->ops->setvalueslocal)(mat,nrow,irow,ncol,icol,y,addv);
1004: }
1005: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1006: return(0);
1007: }
1009: #undef __FUNCT__
1011: /*@
1012: MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
1013: using a local ordering of the nodes a block at a time.
1015: Not Collective
1017: Input Parameters:
1018: + x - the matrix
1019: . nrow, irow - number of rows and their local indices
1020: . ncol, icol - number of columns and their local indices
1021: . y - a logically two-dimensional array of values
1022: - addv - either INSERT_VALUES or ADD_VALUES, where
1023: ADD_VALUES adds values to any existing entries, and
1024: INSERT_VALUES replaces existing entries with new values
1026: Notes:
1027: Before calling MatSetValuesBlockedLocal(), the user must first set the
1028: local-to-global mapping by calling MatSetLocalToGlobalMappingBlock(),
1029: where the mapping MUST be set for matrix blocks, not for matrix elements.
1031: Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1032: options cannot be mixed without intervening calls to the assembly
1033: routines.
1035: These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
1036: MUST be called after all calls to MatSetValuesBlockedLocal() have been completed.
1038: Level: intermediate
1040: Concepts: matrices^putting blocked values in with local numbering
1042: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked()
1043: @*/
1044: int MatSetValuesBlockedLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,PetscScalar *y,InsertMode addv)
1045: {
1046: int ierr,irowm[2048],icolm[2048];
1051: MatPreallocated(mat);
1055: if (mat->insertmode == NOT_SET_VALUES) {
1056: mat->insertmode = addv;
1057: }
1058: #if defined(PETSC_USE_BOPT_g)
1059: else if (mat->insertmode != addv) {
1060: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1061: }
1062: if (!mat->bmapping) {
1063: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with MatSetLocalToGlobalMappingBlock()");
1064: }
1065: if (nrow > 2048 || ncol > 2048) {
1066: SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
1067: }
1068: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1069: #endif
1071: if (mat->assembled) {
1072: mat->was_assembled = PETSC_TRUE;
1073: mat->assembled = PETSC_FALSE;
1074: }
1075: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1076: ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm);
1077: ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm);
1078: (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);
1079: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1080: return(0);
1081: }
1083: /* --------------------------------------------------------*/
1084: #undef __FUNCT__
1086: /*@
1087: MatMult - Computes the matrix-vector product, y = Ax.
1089: Collective on Mat and Vec
1091: Input Parameters:
1092: + mat - the matrix
1093: - x - the vector to be multilplied
1095: Output Parameters:
1096: . y - the result
1098: Notes:
1099: The vectors x and y cannot be the same. I.e., one cannot
1100: call MatMult(A,y,y).
1102: Level: beginner
1104: Concepts: matrix-vector product
1106: .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
1107: @*/
1108: int MatMult(Mat mat,Vec x,Vec y)
1109: {
1115: MatPreallocated(mat);
1119: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1120: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1121: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1122: #ifndef PETSC_HAVE_CONSTRAINTS
1123: if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1124: if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
1125: if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %d %d",mat->m,y->n);
1126: #endif
1128: if (mat->nullsp) {
1129: MatNullSpaceRemove(mat->nullsp,x,&x);
1130: }
1132: PetscLogEventBegin(MAT_Mult,mat,x,y,0);
1133: (*mat->ops->mult)(mat,x,y);
1134: PetscLogEventEnd(MAT_Mult,mat,x,y,0);
1136: if (mat->nullsp) {
1137: MatNullSpaceRemove(mat->nullsp,y,PETSC_NULL);
1138: }
1139: return(0);
1140: }
1142: #undef __FUNCT__
1144: /*@
1145: MatMultTranspose - Computes matrix transpose times a vector.
1147: Collective on Mat and Vec
1149: Input Parameters:
1150: + mat - the matrix
1151: - x - the vector to be multilplied
1153: Output Parameters:
1154: . y - the result
1156: Notes:
1157: The vectors x and y cannot be the same. I.e., one cannot
1158: call MatMultTranspose(A,y,y).
1160: Level: beginner
1162: Concepts: matrix vector product^transpose
1164: .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd()
1165: @*/
1166: int MatMultTranspose(Mat mat,Vec x,Vec y)
1167: {
1169: PetscTruth flg1, flg2;
1174: MatPreallocated(mat);
1178: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1179: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1180: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1181: #ifndef PETSC_HAVE_CONSTRAINTS
1182: if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
1183: if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
1184: #endif
1186: if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP, "Operation not supported");
1187: PetscLogEventBegin(MAT_MultTranspose,mat,x,y,0);
1188: if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply tranpose defined");
1189:
1190: PetscTypeCompare((PetscObject)mat,MATSEQSBAIJ,&flg1);
1191: PetscTypeCompare((PetscObject)mat,MATMPISBAIJ,&flg2);
1192: if (flg1 || flg2) { /* mat is in sbaij format */
1193: (*mat->ops->mult)(mat,x,y);
1194: } else {
1195: (*mat->ops->multtranspose)(mat,x,y);
1196: }
1197: PetscLogEventEnd(MAT_MultTranspose,mat,x,y,0);
1198: return(0);
1199: }
1201: #undef __FUNCT__
1203: /*@
1204: MatMultAdd - Computes v3 = v2 + A * v1.
1206: Collective on Mat and Vec
1208: Input Parameters:
1209: + mat - the matrix
1210: - v1, v2 - the vectors
1212: Output Parameters:
1213: . v3 - the result
1215: Notes:
1216: The vectors v1 and v3 cannot be the same. I.e., one cannot
1217: call MatMultAdd(A,v1,v2,v1).
1219: Level: beginner
1221: Concepts: matrix vector product^addition
1223: .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd()
1224: @*/
1225: int MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1226: {
1232: MatPreallocated(mat);
1237: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1238: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1239: if (mat->N != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->N,v1->N);
1240: if (mat->M != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->M,v2->N);
1241: if (mat->M != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->M,v3->N);
1242: if (mat->m != v3->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: local dim %d %d",mat->m,v3->n);
1243: if (mat->m != v2->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: local dim %d %d",mat->m,v2->n);
1244: if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
1246: PetscLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);
1247: (*mat->ops->multadd)(mat,v1,v2,v3);
1248: PetscLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);
1249: return(0);
1250: }
1252: #undef __FUNCT__
1254: /*@
1255: MatMultTransposeAdd - Computes v3 = v2 + A' * v1.
1257: Collective on Mat and Vec
1259: Input Parameters:
1260: + mat - the matrix
1261: - v1, v2 - the vectors
1263: Output Parameters:
1264: . v3 - the result
1266: Notes:
1267: The vectors v1 and v3 cannot be the same. I.e., one cannot
1268: call MatMultTransposeAdd(A,v1,v2,v1).
1270: Level: beginner
1272: Concepts: matrix vector product^transpose and addition
1274: .seealso: MatMultTranspose(), MatMultAdd(), MatMult()
1275: @*/
1276: int MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1277: {
1283: MatPreallocated(mat);
1288: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1289: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1290: if (!mat->ops->multtransposeadd) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1291: if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
1292: if (mat->M != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->M,v1->N);
1293: if (mat->N != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->N,v2->N);
1294: if (mat->N != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->N,v3->N);
1296: PetscLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);
1297: (*mat->ops->multtransposeadd)(mat,v1,v2,v3);
1298: PetscLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);
1299: return(0);
1300: }
1302: #undef __FUNCT__
1304: /*@
1305: MatMultConstrained - The inner multiplication routine for a
1306: constrained matrix P^T A P.
1308: Collective on Mat and Vec
1310: Input Parameters:
1311: + mat - the matrix
1312: - x - the vector to be multilplied
1314: Output Parameters:
1315: . y - the result
1317: Notes:
1318: The vectors x and y cannot be the same. I.e., one cannot
1319: call MatMult(A,y,y).
1321: Level: beginner
1323: .keywords: matrix, multiply, matrix-vector product, constraint
1324: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
1325: @*/
1326: int MatMultConstrained(Mat mat,Vec x,Vec y)
1327: {
1333: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1334: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1335: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1336: if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1337: if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
1338: if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %d %d",mat->m,y->n);
1340: PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
1341: (*mat->ops->multconstrained)(mat,x,y);
1342: PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
1344: return(0);
1345: }
1347: #undef __FUNCT__
1349: /*@
1350: MatMultTransposeConstrained - The inner multiplication routine for a
1351: constrained matrix P^T A^T P.
1353: Collective on Mat and Vec
1355: Input Parameters:
1356: + mat - the matrix
1357: - x - the vector to be multilplied
1359: Output Parameters:
1360: . y - the result
1362: Notes:
1363: The vectors x and y cannot be the same. I.e., one cannot
1364: call MatMult(A,y,y).
1366: Level: beginner
1368: .keywords: matrix, multiply, matrix-vector product, constraint
1369: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
1370: @*/
1371: int MatMultTransposeConstrained(Mat mat,Vec x,Vec y)
1372: {
1378: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1379: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1380: if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1381: if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1382: if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
1384: PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
1385: (*mat->ops->multtransposeconstrained)(mat,x,y);
1386: PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
1388: return(0);
1389: }
1390: /* ------------------------------------------------------------*/
1391: #undef __FUNCT__
1393: /*@C
1394: MatGetInfo - Returns information about matrix storage (number of
1395: nonzeros, memory, etc.).
1397: Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used
1398: as the flag
1400: Input Parameters:
1401: . mat - the matrix
1403: Output Parameters:
1404: + flag - flag indicating the type of parameters to be returned
1405: (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
1406: MAT_GLOBAL_SUM - sum over all processors)
1407: - info - matrix information context
1409: Notes:
1410: The MatInfo context contains a variety of matrix data, including
1411: number of nonzeros allocated and used, number of mallocs during
1412: matrix assembly, etc. Additional information for factored matrices
1413: is provided (such as the fill ratio, number of mallocs during
1414: factorization, etc.). Much of this info is printed to STDOUT
1415: when using the runtime options
1416: $ -log_info -mat_view_info
1418: Example for C/C++ Users:
1419: See the file ${PETSC_DIR}/include/petscmat.h for a complete list of
1420: data within the MatInfo context. For example,
1421: .vb
1422: MatInfo info;
1423: Mat A;
1424: double mal, nz_a, nz_u;
1426: MatGetInfo(A,MAT_LOCAL,&info);
1427: mal = info.mallocs;
1428: nz_a = info.nz_allocated;
1429: .ve
1431: Example for Fortran Users:
1432: Fortran users should declare info as a double precision
1433: array of dimension MAT_INFO_SIZE, and then extract the parameters
1434: of interest. See the file ${PETSC_DIR}/include/finclude/petscmat.h
1435: a complete list of parameter names.
1436: .vb
1437: double precision info(MAT_INFO_SIZE)
1438: double precision mal, nz_a
1439: Mat A
1440: integer ierr
1442: call MatGetInfo(A,MAT_LOCAL,info,ierr)
1443: mal = info(MAT_INFO_MALLOCS)
1444: nz_a = info(MAT_INFO_NZ_ALLOCATED)
1445: .ve
1447: Level: intermediate
1449: Concepts: matrices^getting information on
1450:
1451: @*/
1452: int MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
1453: {
1459: MatPreallocated(mat);
1461: if (!mat->ops->getinfo) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1462: (*mat->ops->getinfo)(mat,flag,info);
1463: return(0);
1464: }
1466: /* ----------------------------------------------------------*/
1467: #undef __FUNCT__
1469: /*@C
1470: MatILUDTFactor - Performs a drop tolerance ILU factorization.
1472: Collective on Mat
1474: Input Parameters:
1475: + mat - the matrix
1476: . info - information about the factorization to be done
1477: . row - row permutation
1478: - col - column permutation
1480: Output Parameters:
1481: . fact - the factored matrix
1483: Level: developer
1485: Notes:
1486: Most users should employ the simplified SLES interface for linear solvers
1487: instead of working directly with matrix algebra routines such as this.
1488: See, e.g., SLESCreate().
1490: This is currently only supported for the SeqAIJ matrix format using code
1491: from Yousef Saad's SPARSEKIT2 package (translated to C with f2c) and/or
1492: Matlab. SPARSEKIT2 is copyrighted by Yousef Saad with the GNU copyright
1493: and thus can be distributed with PETSc.
1495: Concepts: matrices^ILUDT factorization
1497: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatILUInfo
1498: @*/
1499: int MatILUDTFactor(Mat mat,MatILUInfo *info,IS row,IS col,Mat *fact)
1500: {
1506: MatPreallocated(mat);
1508: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1509: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1510: if (!mat->ops->iludtfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1512: PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
1513: (*mat->ops->iludtfactor)(mat,info,row,col,fact);
1514: PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);
1516: return(0);
1517: }
1519: #undef __FUNCT__
1521: /*@
1522: MatLUFactor - Performs in-place LU factorization of matrix.
1524: Collective on Mat
1526: Input Parameters:
1527: + mat - the matrix
1528: . row - row permutation
1529: . col - column permutation
1530: - info - options for factorization, includes
1531: $ fill - expected fill as ratio of original fill.
1532: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
1533: $ Run with the option -log_info to determine an optimal value to use
1535: Notes:
1536: Most users should employ the simplified SLES interface for linear solvers
1537: instead of working directly with matrix algebra routines such as this.
1538: See, e.g., SLESCreate().
1540: This changes the state of the matrix to a factored matrix; it cannot be used
1541: for example with MatSetValues() unless one first calls MatSetUnfactored().
1543: Level: developer
1545: Concepts: matrices^LU factorization
1547: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(),
1548: MatGetOrdering(), MatSetUnfactored(), MatLUInfo
1550: @*/
1551: int MatLUFactor(Mat mat,IS row,IS col,MatLUInfo *info)
1552: {
1558: MatPreallocated(mat);
1559: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1560: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1561: if (!mat->ops->lufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1563: PetscLogEventBegin(MAT_LUFactor,mat,row,col,0);
1564: (*mat->ops->lufactor)(mat,row,col,info);
1565: PetscLogEventEnd(MAT_LUFactor,mat,row,col,0);
1566: return(0);
1567: }
1569: #undef __FUNCT__
1571: /*@
1572: MatILUFactor - Performs in-place ILU factorization of matrix.
1574: Collective on Mat
1576: Input Parameters:
1577: + mat - the matrix
1578: . row - row permutation
1579: . col - column permutation
1580: - info - structure containing
1581: $ levels - number of levels of fill.
1582: $ expected fill - as ratio of original fill.
1583: $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices
1584: missing diagonal entries)
1586: Notes:
1587: Probably really in-place only when level of fill is zero, otherwise allocates
1588: new space to store factored matrix and deletes previous memory.
1590: Most users should employ the simplified SLES interface for linear solvers
1591: instead of working directly with matrix algebra routines such as this.
1592: See, e.g., SLESCreate().
1594: Level: developer
1596: Concepts: matrices^ILU factorization
1598: .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatILUInfo
1599: @*/
1600: int MatILUFactor(Mat mat,IS row,IS col,MatILUInfo *info)
1601: {
1607: MatPreallocated(mat);
1608: if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
1609: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1610: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1611: if (!mat->ops->ilufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1613: PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
1614: (*mat->ops->ilufactor)(mat,row,col,info);
1615: PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);
1616: return(0);
1617: }
1619: #undef __FUNCT__
1621: /*@
1622: MatLUFactorSymbolic - Performs symbolic LU factorization of matrix.
1623: Call this routine before calling MatLUFactorNumeric().
1625: Collective on Mat
1627: Input Parameters:
1628: + mat - the matrix
1629: . row, col - row and column permutations
1630: - info - options for factorization, includes
1631: $ fill - expected fill as ratio of original fill.
1632: $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
1633: $ Run with the option -log_info to determine an optimal value to use
1635: Output Parameter:
1636: . fact - new matrix that has been symbolically factored
1638: Notes:
1639: See the users manual for additional information about
1640: choosing the fill factor for better efficiency.
1642: Most users should employ the simplified SLES interface for linear solvers
1643: instead of working directly with matrix algebra routines such as this.
1644: See, e.g., SLESCreate().
1646: Level: developer
1648: Concepts: matrices^LU symbolic factorization
1650: .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor(), MatLUInfo
1651: @*/
1652: int MatLUFactorSymbolic(Mat mat,IS row,IS col,MatLUInfo *info,Mat *fact)
1653: {
1659: MatPreallocated(mat);
1663: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1664: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1665: if (!mat->ops->lufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic LU",mat->type_name);
1667: PetscLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);
1668: (*mat->ops->lufactorsymbolic)(mat,row,col,info,fact);
1669: PetscLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);
1670: return(0);
1671: }
1673: #undef __FUNCT__
1675: /*@
1676: MatLUFactorNumeric - Performs numeric LU factorization of a matrix.
1677: Call this routine after first calling MatLUFactorSymbolic().
1679: Collective on Mat
1681: Input Parameters:
1682: + mat - the matrix
1683: - fact - the matrix generated for the factor, from MatLUFactorSymbolic()
1685: Notes:
1686: See MatLUFactor() for in-place factorization. See
1687: MatCholeskyFactorNumeric() for the symmetric, positive definite case.
1689: Most users should employ the simplified SLES interface for linear solvers
1690: instead of working directly with matrix algebra routines such as this.
1691: See, e.g., SLESCreate().
1693: Level: developer
1695: Concepts: matrices^LU numeric factorization
1697: .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
1698: @*/
1699: int MatLUFactorNumeric(Mat mat,Mat *fact)
1700: {
1701: int ierr;
1702: PetscTruth flg;
1707: MatPreallocated(mat);
1710: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1711: if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1712: SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dimensions are different %d should = %d %d should = %d",
1713: mat->M,(*fact)->M,mat->N,(*fact)->N);
1714: }
1715: if (!(*fact)->ops->lufactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1717: PetscLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0);
1718: (*(*fact)->ops->lufactornumeric)(mat,fact);
1719: PetscLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0);
1720: PetscOptionsHasName(PETSC_NULL,"-mat_view_draw",&flg);
1721: if (flg) {
1722: PetscOptionsHasName(PETSC_NULL,"-mat_view_contour",&flg);
1723: if (flg) {
1724: PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);
1725: }
1726: MatView(*fact,PETSC_VIEWER_DRAW_(mat->comm));
1727: PetscViewerFlush(PETSC_VIEWER_DRAW_(mat->comm));
1728: if (flg) {
1729: PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));
1730: }
1731: }
1732: return(0);
1733: }
1735: #undef __FUNCT__
1737: /*@
1738: MatCholeskyFactor - Performs in-place Cholesky factorization of a
1739: symmetric matrix.
1741: Collective on Mat
1743: Input Parameters:
1744: + mat - the matrix
1745: . perm - row and column permutations
1746: - f - expected fill as ratio of original fill
1748: Notes:
1749: See MatLUFactor() for the nonsymmetric case. See also
1750: MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().
1752: Most users should employ the simplified SLES interface for linear solvers
1753: instead of working directly with matrix algebra routines such as this.
1754: See, e.g., SLESCreate().
1756: Level: developer
1758: Concepts: matrices^Cholesky factorization
1760: .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
1761: MatGetOrdering()
1763: @*/
1764: int MatCholeskyFactor(Mat mat,IS perm,PetscReal f)
1765: {
1771: MatPreallocated(mat);
1773: if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
1774: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1775: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1776: if (!mat->ops->choleskyfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1778: PetscLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);
1779: (*mat->ops->choleskyfactor)(mat,perm,f);
1780: PetscLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);
1781: return(0);
1782: }
1784: #undef __FUNCT__
1786: /*@
1787: MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
1788: of a symmetric matrix.
1790: Collective on Mat
1792: Input Parameters:
1793: + mat - the matrix
1794: . perm - row and column permutations
1795: - f - expected fill as ratio of original
1797: Output Parameter:
1798: . fact - the factored matrix
1800: Notes:
1801: See MatLUFactorSymbolic() for the nonsymmetric case. See also
1802: MatCholeskyFactor() and MatCholeskyFactorNumeric().
1804: Most users should employ the simplified SLES interface for linear solvers
1805: instead of working directly with matrix algebra routines such as this.
1806: See, e.g., SLESCreate().
1808: Level: developer
1810: Concepts: matrices^Cholesky symbolic factorization
1812: .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
1813: MatGetOrdering()
1815: @*/
1816: int MatCholeskyFactorSymbolic(Mat mat,IS perm,PetscReal f,Mat *fact)
1817: {
1823: MatPreallocated(mat);
1825: if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
1826: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1827: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1828: if (!mat->ops->choleskyfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1830: PetscLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
1831: (*mat->ops->choleskyfactorsymbolic)(mat,perm,f,fact);
1832: PetscLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
1833: return(0);
1834: }
1836: #undef __FUNCT__
1838: /*@
1839: MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
1840: of a symmetric matrix. Call this routine after first calling
1841: MatCholeskyFactorSymbolic().
1843: Collective on Mat
1845: Input Parameter:
1846: . mat - the initial matrix
1848: Output Parameter:
1849: . fact - the factored matrix
1851: Notes:
1852: Most users should employ the simplified SLES interface for linear solvers
1853: instead of working directly with matrix algebra routines such as this.
1854: See, e.g., SLESCreate().
1856: Level: developer
1858: Concepts: matrices^Cholesky numeric factorization
1860: .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
1861: @*/
1862: int MatCholeskyFactorNumeric(Mat mat,Mat *fact)
1863: {
1869: MatPreallocated(mat);
1871: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1872: if (!(*fact)->ops->choleskyfactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1873: if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1874: SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dim %d should = %d %d should = %d",
1875: mat->M,(*fact)->M,mat->N,(*fact)->N);
1876: }
1878: PetscLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
1879: (*(*fact)->ops->choleskyfactornumeric)(mat,fact);
1880: PetscLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
1881: return(0);
1882: }
1884: /* ----------------------------------------------------------------*/
1885: #undef __FUNCT__
1887: /*@
1888: MatSolve - Solves A x = b, given a factored matrix.
1890: Collective on Mat and Vec
1892: Input Parameters:
1893: + mat - the factored matrix
1894: - b - the right-hand-side vector
1896: Output Parameter:
1897: . x - the result vector
1899: Notes:
1900: The vectors b and x cannot be the same. I.e., one cannot
1901: call MatSolve(A,x,x).
1903: Notes:
1904: Most users should employ the simplified SLES interface for linear solvers
1905: instead of working directly with matrix algebra routines such as this.
1906: See, e.g., SLESCreate().
1908: Level: developer
1910: Concepts: matrices^triangular solves
1912: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd()
1913: @*/
1914: int MatSolve(Mat mat,Vec b,Vec x)
1915: {
1921: MatPreallocated(mat);
1926: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
1927: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
1928: if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1929: if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1930: if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1931: if (mat->M == 0 && mat->N == 0) return(0);
1933: if (!mat->ops->solve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1934: PetscLogEventBegin(MAT_Solve,mat,b,x,0);
1935: (*mat->ops->solve)(mat,b,x);
1936: PetscLogEventEnd(MAT_Solve,mat,b,x,0);
1937: return(0);
1938: }
1940: #undef __FUNCT__
1942: /* @
1943: MatForwardSolve - Solves L x = b, given a factored matrix, A = LU.
1945: Collective on Mat and Vec
1947: Input Parameters:
1948: + mat - the factored matrix
1949: - b - the right-hand-side vector
1951: Output Parameter:
1952: . x - the result vector
1954: Notes:
1955: MatSolve() should be used for most applications, as it performs
1956: a forward solve followed by a backward solve.
1958: The vectors b and x cannot be the same. I.e., one cannot
1959: call MatForwardSolve(A,x,x).
1961: Most users should employ the simplified SLES interface for linear solvers
1962: instead of working directly with matrix algebra routines such as this.
1963: See, e.g., SLESCreate().
1965: Level: developer
1967: Concepts: matrices^forward solves
1969: .seealso: MatSolve(), MatBackwardSolve()
1970: @ */
1971: int MatForwardSolve(Mat mat,Vec b,Vec x)
1972: {
1978: MatPreallocated(mat);
1983: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
1984: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
1985: if (!mat->ops->forwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1986: if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1987: if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1988: if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1990: PetscLogEventBegin(MAT_ForwardSolve,mat,b,x,0);
1991: (*mat->ops->forwardsolve)(mat,b,x);
1992: PetscLogEventEnd(MAT_ForwardSolve,mat,b,x,0);
1993: return(0);
1994: }
1996: #undef __FUNCT__
1998: /* @
1999: MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.
2001: Collective on Mat and Vec
2003: Input Parameters:
2004: + mat - the factored matrix
2005: - b - the right-hand-side vector
2007: Output Parameter:
2008: . x - the result vector
2010: Notes:
2011: MatSolve() should be used for most applications, as it performs
2012: a forward solve followed by a backward solve.
2014: The vectors b and x cannot be the same. I.e., one cannot
2015: call MatBackwardSolve(A,x,x).
2017: Most users should employ the simplified SLES interface for linear solvers
2018: instead of working directly with matrix algebra routines such as this.
2019: See, e.g., SLESCreate().
2021: Level: developer
2023: Concepts: matrices^backward solves
2025: .seealso: MatSolve(), MatForwardSolve()
2026: @ */
2027: int MatBackwardSolve(Mat mat,Vec b,Vec x)
2028: {
2034: MatPreallocated(mat);
2039: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2040: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2041: if (!mat->ops->backwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2042: if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
2043: if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
2044: if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
2046: PetscLogEventBegin(MAT_BackwardSolve,mat,b,x,0);
2047: (*mat->ops->backwardsolve)(mat,b,x);
2048: PetscLogEventEnd(MAT_BackwardSolve,mat,b,x,0);
2049: return(0);
2050: }
2052: #undef __FUNCT__
2054: /*@
2055: MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.
2057: Collective on Mat and Vec
2059: Input Parameters:
2060: + mat - the factored matrix
2061: . b - the right-hand-side vector
2062: - y - the vector to be added to
2064: Output Parameter:
2065: . x - the result vector
2067: Notes:
2068: The vectors b and x cannot be the same. I.e., one cannot
2069: call MatSolveAdd(A,x,y,x).
2071: Most users should employ the simplified SLES interface for linear solvers
2072: instead of working directly with matrix algebra routines such as this.
2073: See, e.g., SLESCreate().
2075: Level: developer
2077: Concepts: matrices^triangular solves
2079: .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd()
2080: @*/
2081: int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
2082: {
2083: PetscScalar one = 1.0;
2084: Vec tmp;
2085: int ierr;
2090: MatPreallocated(mat);
2097: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2098: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2099: if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
2100: if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
2101: if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
2102: if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
2103: if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n);
2105: PetscLogEventBegin(MAT_SolveAdd,mat,b,x,y);
2106: if (mat->ops->solveadd) {
2107: (*mat->ops->solveadd)(mat,b,y,x);
2108: } else {
2109: /* do the solve then the add manually */
2110: if (x != y) {
2111: MatSolve(mat,b,x);
2112: VecAXPY(&one,y,x);
2113: } else {
2114: VecDuplicate(x,&tmp);
2115: PetscLogObjectParent(mat,tmp);
2116: VecCopy(x,tmp);
2117: MatSolve(mat,b,x);
2118: VecAXPY(&one,tmp,x);
2119: VecDestroy(tmp);
2120: }
2121: }
2122: PetscLogEventEnd(MAT_SolveAdd,mat,b,x,y);
2123: return(0);
2124: }
2126: #undef __FUNCT__
2128: /*@
2129: MatSolveTranspose - Solves A' x = b, given a factored matrix.
2131: Collective on Mat and Vec
2133: Input Parameters:
2134: + mat - the factored matrix
2135: - b - the right-hand-side vector
2137: Output Parameter:
2138: . x - the result vector
2140: Notes:
2141: The vectors b and x cannot be the same. I.e., one cannot
2142: call MatSolveTranspose(A,x,x).
2144: Most users should employ the simplified SLES interface for linear solvers
2145: instead of working directly with matrix algebra routines such as this.
2146: See, e.g., SLESCreate().
2148: Level: developer
2150: Concepts: matrices^triangular solves
2152: .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd()
2153: @*/
2154: int MatSolveTranspose(Mat mat,Vec b,Vec x)
2155: {
2161: MatPreallocated(mat);
2166: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2167: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2168: if (!mat->ops->solvetranspose) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s",mat->type_name);
2169: if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
2170: if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);
2172: PetscLogEventBegin(MAT_SolveTranspose,mat,b,x,0);
2173: (*mat->ops->solvetranspose)(mat,b,x);
2174: PetscLogEventEnd(MAT_SolveTranspose,mat,b,x,0);
2175: return(0);
2176: }
2178: #undef __FUNCT__
2180: /*@
2181: MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a
2182: factored matrix.
2184: Collective on Mat and Vec
2186: Input Parameters:
2187: + mat - the factored matrix
2188: . b - the right-hand-side vector
2189: - y - the vector to be added to
2191: Output Parameter:
2192: . x - the result vector
2194: Notes:
2195: The vectors b and x cannot be the same. I.e., one cannot
2196: call MatSolveTransposeAdd(A,x,y,x).
2198: Most users should employ the simplified SLES interface for linear solvers
2199: instead of working directly with matrix algebra routines such as this.
2200: See, e.g., SLESCreate().
2202: Level: developer
2204: Concepts: matrices^triangular solves
2206: .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose()
2207: @*/
2208: int MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x)
2209: {
2210: PetscScalar one = 1.0;
2211: int ierr;
2212: Vec tmp;
2217: MatPreallocated(mat);
2224: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2225: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2226: if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
2227: if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);
2228: if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
2229: if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n);
2231: PetscLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);
2232: if (mat->ops->solvetransposeadd) {
2233: (*mat->ops->solvetransposeadd)(mat,b,y,x);
2234: } else {
2235: /* do the solve then the add manually */
2236: if (x != y) {
2237: MatSolveTranspose(mat,b,x);
2238: VecAXPY(&one,y,x);
2239: } else {
2240: VecDuplicate(x,&tmp);
2241: PetscLogObjectParent(mat,tmp);
2242: VecCopy(x,tmp);
2243: MatSolveTranspose(mat,b,x);
2244: VecAXPY(&one,tmp,x);
2245: VecDestroy(tmp);
2246: }
2247: }
2248: PetscLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);
2249: return(0);
2250: }
2251: /* ----------------------------------------------------------------*/
2253: #undef __FUNCT__
2255: /*@
2256: MatRelax - Computes one relaxation sweep.
2258: Collective on Mat and Vec
2260: Input Parameters:
2261: + mat - the matrix
2262: . b - the right hand side
2263: . omega - the relaxation factor
2264: . flag - flag indicating the type of SOR (see below)
2265: . shift - diagonal shift
2266: - its - the number of iterations
2267: - lits - the number of local iterations
2269: Output Parameters:
2270: . x - the solution (can contain an initial guess)
2272: SOR Flags:
2273: . SOR_FORWARD_SWEEP - forward SOR
2274: . SOR_BACKWARD_SWEEP - backward SOR
2275: . SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
2276: . SOR_LOCAL_FORWARD_SWEEP - local forward SOR
2277: . SOR_LOCAL_BACKWARD_SWEEP - local forward SOR
2278: . SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
2279: . SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies
2280: upper/lower triangular part of matrix to
2281: vector (with omega)
2282: . SOR_ZERO_INITIAL_GUESS - zero initial guess
2284: Notes:
2285: SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
2286: SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings
2287: on each processor.
2289: Application programmers will not generally use MatRelax() directly,
2290: but instead will employ the SLES/PC interface.
2292: Notes for Advanced Users:
2293: The flags are implemented as bitwise inclusive or operations.
2294: For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
2295: to specify a zero initial guess for SSOR.
2297: Most users should employ the simplified SLES interface for linear solvers
2298: instead of working directly with matrix algebra routines such as this.
2299: See, e.g., SLESCreate().
2301: Level: developer
2303: Concepts: matrices^relaxation
2304: Concepts: matrices^SOR
2305: Concepts: matrices^Gauss-Seidel
2307: @*/
2308: int MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,int its,int lits,Vec x)
2309: {
2315: MatPreallocated(mat);
2320: if (!mat->ops->relax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2321: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2322: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2323: if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
2324: if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
2325: if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
2327: PetscLogEventBegin(MAT_Relax,mat,b,x,0);
2328: ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,lits,x);
2329: PetscLogEventEnd(MAT_Relax,mat,b,x,0);
2330: return(0);
2331: }
2333: #undef __FUNCT__
2335: /*
2336: Default matrix copy routine.
2337: */
2338: int MatCopy_Basic(Mat A,Mat B,MatStructure str)
2339: {
2340: int ierr,i,rstart,rend,nz,*cwork;
2341: PetscScalar *vwork;
2344: MatZeroEntries(B);
2345: MatGetOwnershipRange(A,&rstart,&rend);
2346: for (i=rstart; i<rend; i++) {
2347: MatGetRow(A,i,&nz,&cwork,&vwork);
2348: MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);
2349: MatRestoreRow(A,i,&nz,&cwork,&vwork);
2350: }
2351: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2352: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2353: return(0);
2354: }
2356: #undef __FUNCT__
2358: /*@C
2359: MatCopy - Copys a matrix to another matrix.
2361: Collective on Mat
2363: Input Parameters:
2364: + A - the matrix
2365: - str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN
2367: Output Parameter:
2368: . B - where the copy is put
2370: Notes:
2371: If you use SAME_NONZERO_PATTERN then the two matrices had better have the
2372: same nonzero pattern or the routine will crash.
2374: MatCopy() copies the matrix entries of a matrix to another existing
2375: matrix (after first zeroing the second matrix). A related routine is
2376: MatConvert(), which first creates a new matrix and then copies the data.
2378: Level: intermediate
2379:
2380: Concepts: matrices^copying
2382: .seealso: MatConvert(), MatDuplicate()
2384: @*/
2385: int MatCopy(Mat A,Mat B,MatStructure str)
2386: {
2393: MatPreallocated(A);
2395: MatPreallocated(B);
2397: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2398: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2399: if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim (%d,%d) (%d,%d)",A->M,B->M,
2400: A->N,B->N);
2402: PetscLogEventBegin(MAT_Copy,A,B,0,0);
2403: if (A->ops->copy) {
2404: (*A->ops->copy)(A,B,str);
2405: } else { /* generic conversion */
2406: MatCopy_Basic(A,B,str);
2407: }
2408: if (A->mapping) {
2409: if (B->mapping) {ISLocalToGlobalMappingDestroy(B->mapping);B->mapping = 0;}
2410: MatSetLocalToGlobalMapping(B,A->mapping);
2411: }
2412: if (A->bmapping) {
2413: if (B->bmapping) {ISLocalToGlobalMappingDestroy(B->bmapping);B->bmapping = 0;}
2414: MatSetLocalToGlobalMappingBlock(B,A->mapping);
2415: }
2416: PetscLogEventEnd(MAT_Copy,A,B,0,0);
2417: return(0);
2418: }
2420: #include petscsys.h
2421: PetscTruth MatConvertRegisterAllCalled = PETSC_FALSE;
2422: PetscFList MatConvertList = 0;
2424: #undef __FUNCT__
2426: /*@C
2427: MatConvertRegister - Allows one to register a routine that converts a sparse matrix
2428: from one format to another.
2430: Not Collective
2432: Input Parameters:
2433: + type - the type of matrix (defined in include/petscmat.h), for example, MATSEQAIJ.
2434: - Converter - the function that reads the matrix from the binary file.
2436: Level: developer
2438: .seealso: MatConvertRegisterAll(), MatConvert()
2440: @*/
2441: int MatConvertRegister(char *sname,char *path,char *name,int (*function)(Mat,MatType,Mat*))
2442: {
2443: int ierr;
2444: char fullname[PETSC_MAX_PATH_LEN];
2447: PetscFListConcat(path,name,fullname);
2448: PetscFListAdd(&MatConvertList,sname,fullname,(void (*)(void))function);
2449: return(0);
2450: }
2452: #undef __FUNCT__
2454: /*@C
2455: MatConvert - Converts a matrix to another matrix, either of the same
2456: or different type.
2458: Collective on Mat
2460: Input Parameters:
2461: + mat - the matrix
2462: - newtype - new matrix type. Use MATSAME to create a new matrix of the
2463: same type as the original matrix.
2465: Output Parameter:
2466: . M - pointer to place new matrix
2468: Notes:
2469: MatConvert() first creates a new matrix and then copies the data from
2470: the first matrix. A related routine is MatCopy(), which copies the matrix
2471: entries of one matrix to another already existing matrix context.
2473: Level: intermediate
2475: Concepts: matrices^converting between storage formats
2477: .seealso: MatCopy(), MatDuplicate()
2478: @*/
2479: int MatConvert(Mat mat,MatType newtype,Mat *M)
2480: {
2481: int ierr;
2482: PetscTruth sametype,issame,flg;
2483: char convname[256],mtype[256];
2488: MatPreallocated(mat);
2490: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2491: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2493: PetscOptionsGetString(PETSC_NULL,"-matconvert_type",mtype,256,&flg);
2494: if (flg) {
2495: newtype = mtype;
2496: }
2497: PetscLogEventBegin(MAT_Convert,mat,0,0,0);
2498:
2499: PetscTypeCompare((PetscObject)mat,newtype,&sametype);
2500: PetscStrcmp(newtype,"same",&issame);
2501: if ((sametype || issame) && mat->ops->duplicate) {
2502: (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);
2503: } else {
2504: int (*conv)(Mat,MatType,Mat*);
2505: if (!MatConvertRegisterAllCalled) {
2506: MatConvertRegisterAll(PETSC_NULL);
2507: }
2508: PetscFListFind(mat->comm,MatConvertList,newtype,(void(**)(void))&conv);
2509: if (conv) {
2510: (*conv)(mat,newtype,M);
2511: } else {
2512: PetscStrcpy(convname,"MatConvert_");
2513: PetscStrcat(convname,mat->type_name);
2514: PetscStrcat(convname,"_");
2515: PetscStrcat(convname,newtype);
2516: PetscStrcat(convname,"_C");
2517: PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
2518: if (conv) {
2519: (*conv)(mat,newtype,M);
2520: } else {
2521: if (mat->ops->convert) {
2522: (*mat->ops->convert)(mat,newtype,M);
2523: } else {
2524: MatConvert_Basic(mat,newtype,M);
2525: }
2526: }
2527: }
2528: }
2529: PetscLogEventEnd(MAT_Convert,mat,0,0,0);
2530: return(0);
2531: }
2534: #undef __FUNCT__
2536: /*@C
2537: MatDuplicate - Duplicates a matrix including the non-zero structure.
2539: Collective on Mat
2541: Input Parameters:
2542: + mat - the matrix
2543: - op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero
2544: values as well or not
2546: Output Parameter:
2547: . M - pointer to place new matrix
2549: Level: intermediate
2551: Concepts: matrices^duplicating
2553: .seealso: MatCopy(), MatConvert()
2554: @*/
2555: int MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
2556: {
2562: MatPreallocated(mat);
2564: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2565: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2567: *M = 0;
2568: PetscLogEventBegin(MAT_Convert,mat,0,0,0);
2569: if (!mat->ops->duplicate) {
2570: SETERRQ(PETSC_ERR_SUP,"Not written for this matrix type");
2571: }
2572: (*mat->ops->duplicate)(mat,op,M);
2573: if (mat->mapping) {
2574: MatSetLocalToGlobalMapping(*M,mat->mapping);
2575: }
2576: if (mat->bmapping) {
2577: MatSetLocalToGlobalMappingBlock(*M,mat->mapping);
2578: }
2579: PetscLogEventEnd(MAT_Convert,mat,0,0,0);
2580: return(0);
2581: }
2583: #undef __FUNCT__
2585: /*@
2586: MatGetDiagonal - Gets the diagonal of a matrix.
2588: Collective on Mat and Vec
2590: Input Parameters:
2591: + mat - the matrix
2592: - v - the vector for storing the diagonal
2594: Output Parameter:
2595: . v - the diagonal of the matrix
2597: Notes:
2598: For the SeqAIJ matrix format, this routine may also be called
2599: on a LU factored matrix; in that case it routines the reciprocal of
2600: the diagonal entries in U. It returns the entries permuted by the
2601: row and column permutation used during the symbolic factorization.
2603: Level: intermediate
2605: Concepts: matrices^accessing diagonals
2607: .seealso: MatGetRow(), MatGetSubmatrices(), MatGetSubmatrix(), MatGetRowMax()
2608: @*/
2609: int MatGetDiagonal(Mat mat,Vec v)
2610: {
2616: MatPreallocated(mat);
2619: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2620: if (!mat->ops->getdiagonal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2622: (*mat->ops->getdiagonal)(mat,v);
2623: return(0);
2624: }
2626: #undef __FUNCT__
2628: /*@
2629: MatGetRowMax - Gets the maximum value (in absolute value) of each
2630: row of the matrix
2632: Collective on Mat and Vec
2634: Input Parameters:
2635: . mat - the matrix
2637: Output Parameter:
2638: . v - the vector for storing the maximums
2640: Level: intermediate
2642: Concepts: matrices^getting row maximums
2644: .seealso: MatGetDiagonal(), MatGetSubmatrices(), MatGetSubmatrix()
2645: @*/
2646: int MatGetRowMax(Mat mat,Vec v)
2647: {
2653: MatPreallocated(mat);
2656: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2657: if (!mat->ops->getrowmax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2659: (*mat->ops->getrowmax)(mat,v);
2660: return(0);
2661: }
2663: #undef __FUNCT__
2665: /*@C
2666: MatTranspose - Computes an in-place or out-of-place transpose of a matrix.
2668: Collective on Mat
2670: Input Parameter:
2671: . mat - the matrix to transpose
2673: Output Parameters:
2674: . B - the transpose (or pass in PETSC_NULL for an in-place transpose)
2676: Level: intermediate
2678: Concepts: matrices^transposing
2680: .seealso: MatMultTranspose(), MatMultTransposeAdd()
2681: @*/
2682: int MatTranspose(Mat mat,Mat *B)
2683: {
2689: MatPreallocated(mat);
2690: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2691: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2692: if (!mat->ops->transpose) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2694: PetscLogEventBegin(MAT_Transpose,mat,0,0,0);
2695: (*mat->ops->transpose)(mat,B);
2696: PetscLogEventEnd(MAT_Transpose,mat,0,0,0);
2697: return(0);
2698: }
2700: #undef __FUNCT__
2702: /*@C
2703: MatPermute - Creates a new matrix with rows and columns permuted from the
2704: original.
2706: Collective on Mat
2708: Input Parameters:
2709: + mat - the matrix to permute
2710: . row - row permutation, each processor supplies only the permutation for its rows
2711: - col - column permutation, each processor needs the entire column permutation, that is
2712: this is the same size as the total number of columns in the matrix
2714: Output Parameters:
2715: . B - the permuted matrix
2717: Level: advanced
2719: Concepts: matrices^permuting
2721: .seealso: MatGetOrdering()
2722: @*/
2723: int MatPermute(Mat mat,IS row,IS col,Mat *B)
2724: {
2730: MatPreallocated(mat);
2733: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2734: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2735: if (!mat->ops->permute) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2736: (*mat->ops->permute)(mat,row,col,B);
2737: return(0);
2738: }
2740: #undef __FUNCT__
2742: /*@C
2743: MatPermuteSparsify - Creates a new matrix with rows and columns permuted from the
2744: original and sparsified to the prescribed tolerance.
2746: Collective on Mat
2748: Input Parameters:
2749: + A - The matrix to permute
2750: . band - The half-bandwidth of the sparsified matrix, or PETSC_DECIDE
2751: . frac - The half-bandwidth as a fraction of the total size, or 0.0
2752: . tol - The drop tolerance
2753: . rowp - The row permutation
2754: - colp - The column permutation
2756: Output Parameter:
2757: . B - The permuted, sparsified matrix
2759: Level: advanced
2761: Note:
2762: The default behavior (band = PETSC_DECIDE and frac = 0.0) is to
2763: restrict the half-bandwidth of the resulting matrix to 5% of the
2764: total matrix size.
2766: .keywords: matrix, permute, sparsify
2768: .seealso: MatGetOrdering(), MatPermute()
2769: @*/
2770: int MatPermuteSparsify(Mat A, int band, PetscReal frac, PetscReal tol, IS rowp, IS colp, Mat *B)
2771: {
2772: IS irowp, icolp;
2773: int *rows, *cols;
2774: int M, N, locRowStart, locRowEnd;
2775: int nz, newNz;
2776: int *cwork, *cnew;
2777: PetscScalar *vwork, *vnew;
2778: int bw, size;
2779: int row, locRow, newRow, col, newCol;
2780: int ierr;
2786: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
2787: if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2788: if (!A->ops->permutesparsify) {
2789: MatGetSize(A, &M, &N);
2790: MatGetOwnershipRange(A, &locRowStart, &locRowEnd);
2791: ISGetSize(rowp, &size);
2792: if (size != M) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %d for row permutation, should be %d", size, M);
2793: ISGetSize(colp, &size);
2794: if (size != N) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %d for column permutation, should be %d", size, N);
2795: ISInvertPermutation(rowp, 0, &irowp);
2796: ISGetIndices(irowp, &rows);
2797: ISInvertPermutation(colp, 0, &icolp);
2798: ISGetIndices(icolp, &cols);
2799: PetscMalloc(N * sizeof(int), &cnew);
2800: PetscMalloc(N * sizeof(PetscScalar), &vnew);
2802: /* Setup bandwidth to include */
2803: if (band == PETSC_DECIDE) {
2804: if (frac <= 0.0)
2805: bw = (int) (M * 0.05);
2806: else
2807: bw = (int) (M * frac);
2808: } else {
2809: if (band <= 0) SETERRQ(PETSC_ERR_ARG_WRONG, "Bandwidth must be a positive integer");
2810: bw = band;
2811: }
2813: /* Put values into new matrix */
2814: MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, B);
2815: for(row = locRowStart, locRow = 0; row < locRowEnd; row++, locRow++) {
2816: MatGetRow(A, row, &nz, &cwork, &vwork);
2817: newRow = rows[locRow]+locRowStart;
2818: for(col = 0, newNz = 0; col < nz; col++) {
2819: newCol = cols[cwork[col]];
2820: if ((newCol >= newRow - bw) && (newCol < newRow + bw) && (PetscAbsScalar(vwork[col]) >= tol)) {
2821: cnew[newNz] = newCol;
2822: vnew[newNz] = vwork[col];
2823: newNz++;
2824: }
2825: }
2826: MatSetValues(*B, 1, &newRow, newNz, cnew, vnew, INSERT_VALUES);
2827: MatRestoreRow(A, row, &nz, &cwork, &vwork);
2828: }
2829: PetscFree(cnew);
2830: PetscFree(vnew);
2831: MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
2832: MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
2833: ISRestoreIndices(irowp, &rows);
2834: ISRestoreIndices(icolp, &cols);
2835: ISDestroy(irowp);
2836: ISDestroy(icolp);
2837: } else {
2838: (*A->ops->permutesparsify)(A, band, frac, tol, rowp, colp, B);
2839: }
2840: return(0);
2841: }
2843: #undef __FUNCT__
2845: /*@
2846: MatEqual - Compares two matrices.
2848: Collective on Mat
2850: Input Parameters:
2851: + A - the first matrix
2852: - B - the second matrix
2854: Output Parameter:
2855: . flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.
2857: Level: intermediate
2859: Concepts: matrices^equality between
2860: @*/
2861: int MatEqual(Mat A,Mat B,PetscTruth *flg)
2862: {
2869: MatPreallocated(A);
2871: MatPreallocated(B);
2874: if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2875: if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2876: if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim %d %d %d %d",A->M,B->M,A->N,B->N);
2877: if (!A->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",A->type_name);
2878: (*A->ops->equal)(A,B,flg);
2879: return(0);
2880: }
2882: #undef __FUNCT__
2884: /*@
2885: MatDiagonalScale - Scales a matrix on the left and right by diagonal
2886: matrices that are stored as vectors. Either of the two scaling
2887: matrices can be PETSC_NULL.
2889: Collective on Mat
2891: Input Parameters:
2892: + mat - the matrix to be scaled
2893: . l - the left scaling vector (or PETSC_NULL)
2894: - r - the right scaling vector (or PETSC_NULL)
2896: Notes:
2897: MatDiagonalScale() computes A = LAR, where
2898: L = a diagonal matrix, R = a diagonal matrix
2900: Level: intermediate
2902: Concepts: matrices^diagonal scaling
2903: Concepts: diagonal scaling of matrices
2905: .seealso: MatScale()
2906: @*/
2907: int MatDiagonalScale(Mat mat,Vec l,Vec r)
2908: {
2914: MatPreallocated(mat);
2915: if (!mat->ops->diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2918: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2919: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2921: PetscLogEventBegin(MAT_Scale,mat,0,0,0);
2922: (*mat->ops->diagonalscale)(mat,l,r);
2923: PetscLogEventEnd(MAT_Scale,mat,0,0,0);
2924: return(0);
2925: }
2927: #undef __FUNCT__
2929: /*@
2930: MatScale - Scales all elements of a matrix by a given number.
2932: Collective on Mat
2934: Input Parameters:
2935: + mat - the matrix to be scaled
2936: - a - the scaling value
2938: Output Parameter:
2939: . mat - the scaled matrix
2941: Level: intermediate
2943: Concepts: matrices^scaling all entries
2945: .seealso: MatDiagonalScale()
2946: @*/
2947: int MatScale(PetscScalar *a,Mat mat)
2948: {
2954: MatPreallocated(mat);
2956: if (!mat->ops->scale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2957: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2958: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2960: PetscLogEventBegin(MAT_Scale,mat,0,0,0);
2961: (*mat->ops->scale)(a,mat);
2962: PetscLogEventEnd(MAT_Scale,mat,0,0,0);
2963: return(0);
2964: }
2966: #undef __FUNCT__
2968: /*@
2969: MatNorm - Calculates various norms of a matrix.
2971: Collective on Mat
2973: Input Parameters:
2974: + mat - the matrix
2975: - type - the type of norm, NORM_1, NORM_FROBENIUS, NORM_INFINITY
2977: Output Parameters:
2978: . nrm - the resulting norm
2980: Level: intermediate
2982: Concepts: matrices^norm
2983: Concepts: norm^of matrix
2984: @*/
2985: int MatNorm(Mat mat,NormType type,PetscReal *nrm)
2986: {
2992: MatPreallocated(mat);
2995: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2996: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2997: if (!mat->ops->norm) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2998: (*mat->ops->norm)(mat,type,nrm);
2999: return(0);
3000: }
3002: /*
3003: This variable is used to prevent counting of MatAssemblyBegin() that
3004: are called from within a MatAssemblyEnd().
3005: */
3006: static int MatAssemblyEnd_InUse = 0;
3007: #undef __FUNCT__
3009: /*@
3010: MatAssemblyBegin - Begins assembling the matrix. This routine should
3011: be called after completing all calls to MatSetValues().
3013: Collective on Mat
3015: Input Parameters:
3016: + mat - the matrix
3017: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
3018:
3019: Notes:
3020: MatSetValues() generally caches the values. The matrix is ready to
3021: use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
3022: Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
3023: in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
3024: using the matrix.
3026: Level: beginner
3028: Concepts: matrices^assembling
3030: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled()
3031: @*/
3032: int MatAssemblyBegin(Mat mat,MatAssemblyType type)
3033: {
3039: MatPreallocated(mat);
3040: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.nDid you forget to call MatSetUnfactored()?");
3041: if (mat->assembled) {
3042: mat->was_assembled = PETSC_TRUE;
3043: mat->assembled = PETSC_FALSE;
3044: }
3045: if (!MatAssemblyEnd_InUse) {
3046: PetscLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);
3047: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
3048: PetscLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);
3049: } else {
3050: if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
3051: }
3052: return(0);
3053: }
3055: #undef __FUNCT__
3057: /*@
3058: MatAssembled - Indicates if a matrix has been assembled and is ready for
3059: use; for example, in matrix-vector product.
3061: Collective on Mat
3063: Input Parameter:
3064: . mat - the matrix
3066: Output Parameter:
3067: . assembled - PETSC_TRUE or PETSC_FALSE
3069: Level: advanced
3071: Concepts: matrices^assembled?
3073: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin()
3074: @*/
3075: int MatAssembled(Mat mat,PetscTruth *assembled)
3076: {
3080: MatPreallocated(mat);
3081: *assembled = mat->assembled;
3082: return(0);
3083: }
3085: #undef __FUNCT__
3087: /*
3088: Processes command line options to determine if/how a matrix
3089: is to be viewed. Called by MatAssemblyEnd() and MatLoad().
3090: */
3091: int MatView_Private(Mat mat)
3092: {
3093: int ierr;
3094: PetscTruth flg;
3095: static PetscTruth incall = PETSC_FALSE;
3098: if (incall) return(0);
3099: incall = PETSC_TRUE;
3100: PetscOptionsBegin(mat->comm,mat->prefix,"Matrix Options","Mat");
3101: PetscOptionsName("-mat_view_info","Information on matrix size","MatView",&flg);
3102: if (flg) {
3103: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO);
3104: MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
3105: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));
3106: }
3107: PetscOptionsName("-mat_view_info_detailed","Nonzeros in the matrix","MatView",&flg);
3108: if (flg) {
3109: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO_DETAIL);
3110: MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
3111: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));
3112: }
3113: PetscOptionsName("-mat_view","Print matrix to stdout","MatView",&flg);
3114: if (flg) {
3115: MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
3116: }
3117: PetscOptionsName("-mat_view_matlab","Print matrix to stdout in a format Matlab can read","MatView",&flg);
3118: if (flg) {
3119: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_MATLAB);
3120: MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
3121: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));
3122: }
3123: PetscOptionsName("-mat_view_socket","Send matrix to socket (can be read from matlab)","MatView",&flg);
3124: if (flg) {
3125: MatView(mat,PETSC_VIEWER_SOCKET_(mat->comm));
3126: PetscViewerFlush(PETSC_VIEWER_SOCKET_(mat->comm));
3127: }
3128: PetscOptionsName("-mat_view_binary","Save matrix to file in binary format","MatView",&flg);
3129: if (flg) {
3130: MatView(mat,PETSC_VIEWER_BINARY_(mat->comm));
3131: PetscViewerFlush(PETSC_VIEWER_BINARY_(mat->comm));
3132: }
3133: PetscOptionsEnd();
3134: /* cannot have inside PetscOptionsBegin() because uses PetscOptionsBegin() */
3135: PetscOptionsHasName(mat->prefix,"-mat_view_draw",&flg);
3136: if (flg) {
3137: PetscOptionsHasName(mat->prefix,"-mat_view_contour",&flg);
3138: if (flg) {
3139: PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);
3140: }
3141: MatView(mat,PETSC_VIEWER_DRAW_(mat->comm));
3142: PetscViewerFlush(PETSC_VIEWER_DRAW_(mat->comm));
3143: if (flg) {
3144: PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));
3145: }
3146: }
3147: incall = PETSC_FALSE;
3148: return(0);
3149: }
3151: #undef __FUNCT__
3153: /*@
3154: MatAssemblyEnd - Completes assembling the matrix. This routine should
3155: be called after MatAssemblyBegin().
3157: Collective on Mat
3159: Input Parameters:
3160: + mat - the matrix
3161: - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
3163: Options Database Keys:
3164: + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
3165: . -mat_view_info_detailed - Prints more detailed info
3166: . -mat_view - Prints matrix in ASCII format
3167: . -mat_view_matlab - Prints matrix in Matlab format
3168: . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
3169: . -display <name> - Sets display name (default is host)
3170: . -draw_pause <sec> - Sets number of seconds to pause after display
3171: . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
3172: . -viewer_socket_machine <machine>
3173: . -viewer_socket_port <port>
3174: . -mat_view_binary - save matrix to file in binary format
3175: - -viewer_binary_filename <name>
3177: Notes:
3178: MatSetValues() generally caches the values. The matrix is ready to
3179: use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
3180: Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
3181: in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
3182: using the matrix.
3184: Level: beginner
3186: .seealso: MatAssemblyBegin(), MatSetValues(), PetscDrawOpenX(), MatView(), MatAssembled(), PetscViewerSocketOpen()
3187: @*/
3188: int MatAssemblyEnd(Mat mat,MatAssemblyType type)
3189: {
3190: int ierr;
3191: static int inassm = 0;
3192: PetscTruth flg;
3197: MatPreallocated(mat);
3199: inassm++;
3200: MatAssemblyEnd_InUse++;
3201: if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
3202: PetscLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);
3203: if (mat->ops->assemblyend) {
3204: (*mat->ops->assemblyend)(mat,type);
3205: }
3206: PetscLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);
3207: } else {
3208: if (mat->ops->assemblyend) {
3209: (*mat->ops->assemblyend)(mat,type);
3210: }
3211: }
3213: /* Flush assembly is not a true assembly */
3214: if (type != MAT_FLUSH_ASSEMBLY) {
3215: mat->assembled = PETSC_TRUE; mat->num_ass++;
3216: }
3217: mat->insertmode = NOT_SET_VALUES;
3218: MatAssemblyEnd_InUse--;
3220: if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) {
3221: MatView_Private(mat);
3222: }
3223: inassm--;
3224: PetscOptionsHasName(mat->prefix,"-help",&flg);
3225: if (flg) {
3226: MatPrintHelp(mat);
3227: }
3228: return(0);
3229: }
3232: #undef __FUNCT__
3234: /*@
3235: MatCompress - Tries to store the matrix in as little space as
3236: possible. May fail if memory is already fully used, since it
3237: tries to allocate new space.
3239: Collective on Mat
3241: Input Parameters:
3242: . mat - the matrix
3244: Level: advanced
3246: @*/
3247: int MatCompress(Mat mat)
3248: {
3254: MatPreallocated(mat);
3255: if (mat->ops->compress) {(*mat->ops->compress)(mat);}
3256: return(0);
3257: }
3259: #undef __FUNCT__
3261: /*@
3262: MatSetOption - Sets a parameter option for a matrix. Some options
3263: may be specific to certain storage formats. Some options
3264: determine how values will be inserted (or added). Sorted,
3265: row-oriented input will generally assemble the fastest. The default
3266: is row-oriented, nonsorted input.
3268: Collective on Mat
3270: Input Parameters:
3271: + mat - the matrix
3272: - option - the option, one of those listed below (and possibly others),
3273: e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR
3275: Options Describing Matrix Structure:
3276: + MAT_SYMMETRIC - symmetric in terms of both structure and value
3277: - MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure
3279: Options For Use with MatSetValues():
3280: Insert a logically dense subblock, which can be
3281: + MAT_ROW_ORIENTED - row-oriented (default)
3282: . MAT_COLUMN_ORIENTED - column-oriented
3283: . MAT_ROWS_SORTED - sorted by row
3284: . MAT_ROWS_UNSORTED - not sorted by row (default)
3285: . MAT_COLUMNS_SORTED - sorted by column
3286: - MAT_COLUMNS_UNSORTED - not sorted by column (default)
3288: Not these options reflect the data you pass in with MatSetValues(); it has
3289: nothing to do with how the data is stored internally in the matrix
3290: data structure.
3292: When (re)assembling a matrix, we can restrict the input for
3293: efficiency/debugging purposes. These options include
3294: + MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be
3295: allowed if they generate a new nonzero
3296: . MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed
3297: . MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if
3298: they generate a nonzero in a new diagonal (for block diagonal format only)
3299: . MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only)
3300: . MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries
3301: . MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry
3302: - MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly
3304: Notes:
3305: Some options are relevant only for particular matrix types and
3306: are thus ignored by others. Other options are not supported by
3307: certain matrix types and will generate an error message if set.
3309: If using a Fortran 77 module to compute a matrix, one may need to
3310: use the column-oriented option (or convert to the row-oriented
3311: format).
3313: MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion
3314: that would generate a new entry in the nonzero structure is instead
3315: ignored. Thus, if memory has not alredy been allocated for this particular
3316: data, then the insertion is ignored. For dense matrices, in which
3317: the entire array is allocated, no entries are ever ignored.
3318: Set after the first MatAssemblyEnd()
3320: MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion
3321: that would generate a new entry in the nonzero structure instead produces
3322: an error. (Currently supported for AIJ and BAIJ formats only.)
3323: This is a useful flag when using SAME_NONZERO_PATTERN in calling
3324: SLESSetOperators() to ensure that the nonzero pattern truely does
3325: remain unchanged. Set after the first MatAssemblyEnd()
3327: MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion
3328: that would generate a new entry that has not been preallocated will
3329: instead produce an error. (Currently supported for AIJ and BAIJ formats
3330: only.) This is a useful flag when debugging matrix memory preallocation.
3332: MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for
3333: other processors should be dropped, rather than stashed.
3334: This is useful if you know that the "owning" processor is also
3335: always generating the correct matrix entries, so that PETSc need
3336: not transfer duplicate entries generated on another processor.
3337:
3338: MAT_USE_HASH_TABLE indicates that a hash table be used to improve the
3339: searches during matrix assembly. When this flag is set, the hash table
3340: is created during the first Matrix Assembly. This hash table is
3341: used the next time through, during MatSetVaules()/MatSetVaulesBlocked()
3342: to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag
3343: should be used with MAT_USE_HASH_TABLE flag. This option is currently
3344: supported by MATMPIBAIJ format only.
3346: MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries
3347: are kept in the nonzero structure
3349: MAT_IGNORE_ZERO_ENTRIES - when using ADD_VALUES for AIJ matrices this will stop
3350: zero values from creating a zero location in the matrix
3352: MAT_USE_INODES - indicates using inode version of the code - works with AIJ and
3353: ROWBS matrix types
3355: MAT_DO_NOT_USE_INODES - indicates not using inode version of the code - works
3356: with AIJ and ROWBS matrix types
3358: Level: intermediate
3360: Concepts: matrices^setting options
3362: @*/
3363: int MatSetOption(Mat mat,MatOption op)
3364: {
3370: MatPreallocated(mat);
3371: switch (op) {
3372: case MAT_SYMMETRIC:
3373: mat->symmetric = PETSC_TRUE;
3374: mat->structurally_symmetric = PETSC_TRUE;
3375: break;
3376: case MAT_STRUCTURALLY_SYMMETRIC:
3377: mat->structurally_symmetric = PETSC_TRUE;
3378: break;
3379: default:
3380: if (mat->ops->setoption) {(*mat->ops->setoption)(mat,op);}
3381: break;
3382: }
3383: return(0);
3384: }
3386: #undef __FUNCT__
3388: /*@
3389: MatZeroEntries - Zeros all entries of a matrix. For sparse matrices
3390: this routine retains the old nonzero structure.
3392: Collective on Mat
3394: Input Parameters:
3395: . mat - the matrix
3397: Level: intermediate
3399: Concepts: matrices^zeroing
3401: .seealso: MatZeroRows()
3402: @*/
3403: int MatZeroEntries(Mat mat)
3404: {
3410: MatPreallocated(mat);
3411: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3412: if (!mat->ops->zeroentries) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3414: PetscLogEventBegin(MAT_ZeroEntries,mat,0,0,0);
3415: (*mat->ops->zeroentries)(mat);
3416: PetscLogEventEnd(MAT_ZeroEntries,mat,0,0,0);
3417: return(0);
3418: }
3420: #undef __FUNCT__
3422: /*@C
3423: MatZeroRows - Zeros all entries (except possibly the main diagonal)
3424: of a set of rows of a matrix.
3426: Collective on Mat
3428: Input Parameters:
3429: + mat - the matrix
3430: . is - index set of rows to remove
3431: - diag - pointer to value put in all diagonals of eliminated rows.
3432: Note that diag is not a pointer to an array, but merely a
3433: pointer to a single value.
3435: Notes:
3436: For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
3437: but does not release memory. For the dense and block diagonal
3438: formats this does not alter the nonzero structure.
3440: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
3441: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
3442: merely zeroed.
3444: The user can set a value in the diagonal entry (or for the AIJ and
3445: row formats can optionally remove the main diagonal entry from the
3446: nonzero structure as well, by passing a null pointer (PETSC_NULL
3447: in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
3449: For the parallel case, all processes that share the matrix (i.e.,
3450: those in the communicator used for matrix creation) MUST call this
3451: routine, regardless of whether any rows being zeroed are owned by
3452: them.
3454: For the SBAIJ matrix (since only the upper triangular half of the matrix
3455: is stored) the effect of this call is to also zero the corresponding
3456: column.
3457:
3458: Level: intermediate
3460: Concepts: matrices^zeroing rows
3462: .seealso: MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
3463: @*/
3464: int MatZeroRows(Mat mat,IS is,PetscScalar *diag)
3465: {
3471: MatPreallocated(mat);
3474: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3475: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3476: if (!mat->ops->zerorows) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3478: (*mat->ops->zerorows)(mat,is,diag);
3479: MatView_Private(mat);
3480: return(0);
3481: }
3483: #undef __FUNCT__
3485: /*@C
3486: MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
3487: of a set of rows of a matrix; using local numbering of rows.
3489: Collective on Mat
3491: Input Parameters:
3492: + mat - the matrix
3493: . is - index set of rows to remove
3494: - diag - pointer to value put in all diagonals of eliminated rows.
3495: Note that diag is not a pointer to an array, but merely a
3496: pointer to a single value.
3498: Notes:
3499: Before calling MatZeroRowsLocal(), the user must first set the
3500: local-to-global mapping by calling MatSetLocalToGlobalMapping().
3502: For the AIJ matrix formats this removes the old nonzero structure,
3503: but does not release memory. For the dense and block diagonal
3504: formats this does not alter the nonzero structure.
3506: If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
3507: of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
3508: merely zeroed.
3510: The user can set a value in the diagonal entry (or for the AIJ and
3511: row formats can optionally remove the main diagonal entry from the
3512: nonzero structure as well, by passing a null pointer (PETSC_NULL
3513: in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
3515: Level: intermediate
3517: Concepts: matrices^zeroing
3519: .seealso: MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
3520: @*/
3521: int MatZeroRowsLocal(Mat mat,IS is,PetscScalar *diag)
3522: {
3524: IS newis;
3529: MatPreallocated(mat);
3532: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3533: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3535: if (mat->ops->zerorowslocal) {
3536: (*mat->ops->zerorowslocal)(mat,is,diag);
3537: } else {
3538: if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Need to provide local to global mapping to matrix first");
3539: ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);
3540: (*mat->ops->zerorows)(mat,newis,diag);
3541: ISDestroy(newis);
3542: }
3543: return(0);
3544: }
3546: #undef __FUNCT__
3548: /*@
3549: MatGetSize - Returns the numbers of rows and columns in a matrix.
3551: Not Collective
3553: Input Parameter:
3554: . mat - the matrix
3556: Output Parameters:
3557: + m - the number of global rows
3558: - n - the number of global columns
3560: Level: beginner
3562: Concepts: matrices^size
3564: .seealso: MatGetLocalSize()
3565: @*/
3566: int MatGetSize(Mat mat,int *m,int* n)
3567: {
3570: if (m) *m = mat->M;
3571: if (n) *n = mat->N;
3572: return(0);
3573: }
3575: #undef __FUNCT__
3577: /*@
3578: MatGetLocalSize - Returns the number of rows and columns in a matrix
3579: stored locally. This information may be implementation dependent, so
3580: use with care.
3582: Not Collective
3584: Input Parameters:
3585: . mat - the matrix
3587: Output Parameters:
3588: + m - the number of local rows
3589: - n - the number of local columns
3591: Level: beginner
3593: Concepts: matrices^local size
3595: .seealso: MatGetSize()
3596: @*/
3597: int MatGetLocalSize(Mat mat,int *m,int* n)
3598: {
3601: if (m) *m = mat->m;
3602: if (n) *n = mat->n;
3603: return(0);
3604: }
3606: #undef __FUNCT__
3608: /*@
3609: MatGetOwnershipRange - Returns the range of matrix rows owned by
3610: this processor, assuming that the matrix is laid out with the first
3611: n1 rows on the first processor, the next n2 rows on the second, etc.
3612: For certain parallel layouts this range may not be well defined.
3614: Not Collective
3616: Input Parameters:
3617: . mat - the matrix
3619: Output Parameters:
3620: + m - the global index of the first local row
3621: - n - one more than the global index of the last local row
3623: Level: beginner
3625: Concepts: matrices^row ownership
3626: @*/
3627: int MatGetOwnershipRange(Mat mat,int *m,int* n)
3628: {
3634: MatPreallocated(mat);
3637: PetscMapGetLocalRange(mat->rmap,m,n);
3638: return(0);
3639: }
3641: #undef __FUNCT__
3643: /*@
3644: MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
3645: Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric()
3646: to complete the factorization.
3648: Collective on Mat
3650: Input Parameters:
3651: + mat - the matrix
3652: . row - row permutation
3653: . column - column permutation
3654: - info - structure containing
3655: $ levels - number of levels of fill.
3656: $ expected fill - as ratio of original fill.
3657: $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices
3658: missing diagonal entries)
3660: Output Parameters:
3661: . fact - new matrix that has been symbolically factored
3663: Notes:
3664: See the users manual for additional information about
3665: choosing the fill factor for better efficiency.
3667: Most users should employ the simplified SLES interface for linear solvers
3668: instead of working directly with matrix algebra routines such as this.
3669: See, e.g., SLESCreate().
3671: Level: developer
3673: Concepts: matrices^symbolic LU factorization
3674: Concepts: matrices^factorization
3675: Concepts: LU^symbolic factorization
3677: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
3678: MatGetOrdering(), MatILUInfo
3680: @*/
3681: int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatILUInfo *info,Mat *fact)
3682: {
3688: MatPreallocated(mat);
3692: if (info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels of fill negative %d",(int)info->levels);
3693: if (info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",info->fill);
3694: if (!mat->ops->ilufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic ILU",mat->type_name);
3695: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3696: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3698: PetscLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);
3699: (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);
3700: PetscLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);
3701: return(0);
3702: }
3704: #undef __FUNCT__
3706: /*@
3707: MatICCFactorSymbolic - Performs symbolic incomplete
3708: Cholesky factorization for a symmetric matrix. Use
3709: MatCholeskyFactorNumeric() to complete the factorization.
3711: Collective on Mat
3713: Input Parameters:
3714: + mat - the matrix
3715: . perm - row and column permutation
3716: . fill - levels of fill
3717: - f - expected fill as ratio of original fill
3719: Output Parameter:
3720: . fact - the factored matrix
3722: Notes:
3723: Currently only no-fill factorization is supported.
3725: Most users should employ the simplified SLES interface for linear solvers
3726: instead of working directly with matrix algebra routines such as this.
3727: See, e.g., SLESCreate().
3729: Level: developer
3731: Concepts: matrices^symbolic incomplete Cholesky factorization
3732: Concepts: matrices^factorization
3733: Concepts: Cholsky^symbolic factorization
3735: .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor()
3736: @*/
3737: int MatICCFactorSymbolic(Mat mat,IS perm,PetscReal f,int fill,Mat *fact)
3738: {
3744: MatPreallocated(mat);
3747: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3748: if (fill < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Fill negative %d",fill);
3749: if (f < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",f);
3750: if (!mat->ops->iccfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic ICC",mat->type_name);
3751: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3753: PetscLogEventBegin(MAT_ICCFactorSymbolic,mat,perm,0,0);
3754: (*mat->ops->iccfactorsymbolic)(mat,perm,f,fill,fact);
3755: PetscLogEventEnd(MAT_ICCFactorSymbolic,mat,perm,0,0);
3756: return(0);
3757: }
3759: #undef __FUNCT__
3761: /*@C
3762: MatGetArray - Returns a pointer to the element values in the matrix.
3763: The result of this routine is dependent on the underlying matrix data
3764: structure, and may not even work for certain matrix types. You MUST
3765: call MatRestoreArray() when you no longer need to access the array.
3767: Not Collective
3769: Input Parameter:
3770: . mat - the matrix
3772: Output Parameter:
3773: . v - the location of the values
3776: Fortran Note:
3777: This routine is used differently from Fortran, e.g.,
3778: .vb
3779: Mat mat
3780: PetscScalar mat_array(1)
3781: PetscOffset i_mat
3782: int ierr
3783: call MatGetArray(mat,mat_array,i_mat,ierr)
3785: C Access first local entry in matrix; note that array is
3786: C treated as one dimensional
3787: value = mat_array(i_mat + 1)
3789: [... other code ...]
3790: call MatRestoreArray(mat,mat_array,i_mat,ierr)
3791: .ve
3793: See the Fortran chapter of the users manual and
3794: petsc/src/mat/examples/tests for details.
3796: Level: advanced
3798: Concepts: matrices^access array
3800: .seealso: MatRestoreArray(), MatGetArrayF90()
3801: @*/
3802: int MatGetArray(Mat mat,PetscScalar **v)
3803: {
3809: MatPreallocated(mat);
3811: if (!mat->ops->getarray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3812: (*mat->ops->getarray)(mat,v);
3813: return(0);
3814: }
3816: #undef __FUNCT__
3818: /*@C
3819: MatRestoreArray - Restores the matrix after MatGetArray() has been called.
3821: Not Collective
3823: Input Parameter:
3824: + mat - the matrix
3825: - v - the location of the values
3827: Fortran Note:
3828: This routine is used differently from Fortran, e.g.,
3829: .vb
3830: Mat mat
3831: PetscScalar mat_array(1)
3832: PetscOffset i_mat
3833: int ierr
3834: call MatGetArray(mat,mat_array,i_mat,ierr)
3836: C Access first local entry in matrix; note that array is
3837: C treated as one dimensional
3838: value = mat_array(i_mat + 1)
3840: [... other code ...]
3841: call MatRestoreArray(mat,mat_array,i_mat,ierr)
3842: .ve
3844: See the Fortran chapter of the users manual and
3845: petsc/src/mat/examples/tests for details
3847: Level: advanced
3849: .seealso: MatGetArray(), MatRestoreArrayF90()
3850: @*/
3851: int MatRestoreArray(Mat mat,PetscScalar **v)
3852: {
3858: MatPreallocated(mat);
3860: if (!mat->ops->restorearray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3861: (*mat->ops->restorearray)(mat,v);
3862: return(0);
3863: }
3865: #undef __FUNCT__
3867: /*@C
3868: MatGetSubMatrices - Extracts several submatrices from a matrix. If submat
3869: points to an array of valid matrices, they may be reused to store the new
3870: submatrices.
3872: Collective on Mat
3874: Input Parameters:
3875: + mat - the matrix
3876: . n - the number of submatrixes to be extracted (on this processor, may be zero)
3877: . irow, icol - index sets of rows and columns to extract
3878: - scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3880: Output Parameter:
3881: . submat - the array of submatrices
3883: Notes:
3884: MatGetSubMatrices() can extract only sequential submatrices
3885: (from both sequential and parallel matrices). Use MatGetSubMatrix()
3886: to extract a parallel submatrix.
3888: When extracting submatrices from a parallel matrix, each processor can
3889: form a different submatrix by setting the rows and columns of its
3890: individual index sets according to the local submatrix desired.
3892: When finished using the submatrices, the user should destroy
3893: them with MatDestroyMatrices().
3895: MAT_REUSE_MATRIX can only be used when the nonzero structure of the
3896: original matrix has not changed from that last call to MatGetSubMatrices().
3898: This routine creates the matrices in submat; you should NOT create them before
3899: calling it. It also allocates the array of matrix pointers submat.
3901: Fortran Note:
3902: The Fortran interface is slightly different from that given below; it
3903: requires one to pass in as submat a Mat (integer) array of size at least m.
3905: Level: advanced
3907: Concepts: matrices^accessing submatrices
3908: Concepts: submatrices
3910: .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal()
3911: @*/
3912: int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatReuse scall,Mat **submat)
3913: {
3914: int ierr;
3919: MatPreallocated(mat);
3920: if (!mat->ops->getsubmatrices) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3921: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3923: PetscLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);
3924: (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);
3925: PetscLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);
3926: return(0);
3927: }
3929: #undef __FUNCT__
3931: /*@C
3932: MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices().
3934: Collective on Mat
3936: Input Parameters:
3937: + n - the number of local matrices
3938: - mat - the matrices
3940: Level: advanced
3942: Notes: Frees not only the matrices, but also the array that contains the matrices
3944: .seealso: MatGetSubMatrices()
3945: @*/
3946: int MatDestroyMatrices(int n,Mat **mat)
3947: {
3948: int ierr,i;
3951: if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of matrices %d",n);
3953: for (i=0; i<n; i++) {
3954: MatDestroy((*mat)[i]);
3955: }
3956: /* memory is allocated even if n = 0 */
3957: PetscFree(*mat);
3958: return(0);
3959: }
3961: #undef __FUNCT__
3963: /*@
3964: MatIncreaseOverlap - Given a set of submatrices indicated by index sets,
3965: replaces the index sets by larger ones that represent submatrices with
3966: additional overlap.
3968: Collective on Mat
3970: Input Parameters:
3971: + mat - the matrix
3972: . n - the number of index sets
3973: . is - the array of pointers to index sets
3974: - ov - the additional overlap requested
3976: Level: developer
3978: Concepts: overlap
3979: Concepts: ASM^computing overlap
3981: .seealso: MatGetSubMatrices()
3982: @*/
3983: int MatIncreaseOverlap(Mat mat,int n,IS *is,int ov)
3984: {
3990: MatPreallocated(mat);
3991: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3992: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3994: if (!ov) return(0);
3995: if (!mat->ops->increaseoverlap) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3996: PetscLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);
3997: (*mat->ops->increaseoverlap)(mat,n,is,ov);
3998: PetscLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);
3999: return(0);
4000: }
4002: #undef __FUNCT__
4004: /*@
4005: MatPrintHelp - Prints all the options for the matrix.
4007: Collective on Mat
4009: Input Parameter:
4010: . mat - the matrix
4012: Options Database Keys:
4013: + -help - Prints matrix options
4014: - -h - Prints matrix options
4016: Level: developer
4018: .seealso: MatCreate(), MatCreateXXX()
4019: @*/
4020: int MatPrintHelp(Mat mat)
4021: {
4022: static PetscTruth called = PETSC_FALSE;
4023: int ierr;
4028: MatPreallocated(mat);
4030: if (!called) {
4031: if (mat->ops->printhelp) {
4032: (*mat->ops->printhelp)(mat);
4033: }
4034: called = PETSC_TRUE;
4035: }
4036: return(0);
4037: }
4039: #undef __FUNCT__
4041: /*@
4042: MatGetBlockSize - Returns the matrix block size; useful especially for the
4043: block row and block diagonal formats.
4044:
4045: Not Collective
4047: Input Parameter:
4048: . mat - the matrix
4050: Output Parameter:
4051: . bs - block size
4053: Notes:
4054: Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG.
4055: Block row formats are MATSEQBAIJ, MATMPIBAIJ
4057: Level: intermediate
4059: Concepts: matrices^block size
4061: .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag()
4062: @*/
4063: int MatGetBlockSize(Mat mat,int *bs)
4064: {
4070: MatPreallocated(mat);
4072: if (!mat->ops->getblocksize) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4073: (*mat->ops->getblocksize)(mat,bs);
4074: return(0);
4075: }
4077: #undef __FUNCT__
4079: /*@C
4080: MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices.
4082: Collective on Mat
4084: Input Parameters:
4085: + mat - the matrix
4086: . shift - 0 or 1 indicating we want the indices starting at 0 or 1
4087: - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
4088: symmetrized
4090: Output Parameters:
4091: + n - number of rows in the (possibly compressed) matrix
4092: . ia - the row pointers
4093: . ja - the column indices
4094: - done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
4096: Level: developer
4098: .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
4099: @*/
4100: int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
4101: {
4107: MatPreallocated(mat);
4111: if (!mat->ops->getrowij) *done = PETSC_FALSE;
4112: else {
4113: *done = PETSC_TRUE;
4114: ierr = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);
4115: }
4116: return(0);
4117: }
4119: #undef __FUNCT__
4121: /*@C
4122: MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices.
4124: Collective on Mat
4126: Input Parameters:
4127: + mat - the matrix
4128: . shift - 1 or zero indicating we want the indices starting at 0 or 1
4129: - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
4130: symmetrized
4132: Output Parameters:
4133: + n - number of columns in the (possibly compressed) matrix
4134: . ia - the column pointers
4135: . ja - the row indices
4136: - done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
4138: Level: developer
4140: .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
4141: @*/
4142: int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
4143: {
4149: MatPreallocated(mat);
4154: if (!mat->ops->getcolumnij) *done = PETSC_FALSE;
4155: else {
4156: *done = PETSC_TRUE;
4157: ierr = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);
4158: }
4159: return(0);
4160: }
4162: #undef __FUNCT__
4164: /*@C
4165: MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with
4166: MatGetRowIJ().
4168: Collective on Mat
4170: Input Parameters:
4171: + mat - the matrix
4172: . shift - 1 or zero indicating we want the indices starting at 0 or 1
4173: - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
4174: symmetrized
4176: Output Parameters:
4177: + n - size of (possibly compressed) matrix
4178: . ia - the row pointers
4179: . ja - the column indices
4180: - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
4182: Level: developer
4184: .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
4185: @*/
4186: int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
4187: {
4193: MatPreallocated(mat);
4198: if (!mat->ops->restorerowij) *done = PETSC_FALSE;
4199: else {
4200: *done = PETSC_TRUE;
4201: ierr = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);
4202: }
4203: return(0);
4204: }
4206: #undef __FUNCT__
4208: /*@C
4209: MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with
4210: MatGetColumnIJ().
4212: Collective on Mat
4214: Input Parameters:
4215: + mat - the matrix
4216: . shift - 1 or zero indicating we want the indices starting at 0 or 1
4217: - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
4218: symmetrized
4220: Output Parameters:
4221: + n - size of (possibly compressed) matrix
4222: . ia - the column pointers
4223: . ja - the row indices
4224: - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
4226: Level: developer
4228: .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
4229: @*/
4230: int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
4231: {
4237: MatPreallocated(mat);
4242: if (!mat->ops->restorecolumnij) *done = PETSC_FALSE;
4243: else {
4244: *done = PETSC_TRUE;
4245: ierr = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);
4246: }
4247: return(0);
4248: }
4250: #undef __FUNCT__
4252: /*@C
4253: MatColoringPatch -Used inside matrix coloring routines that
4254: use MatGetRowIJ() and/or MatGetColumnIJ().
4256: Collective on Mat
4258: Input Parameters:
4259: + mat - the matrix
4260: . n - number of colors
4261: - colorarray - array indicating color for each column
4263: Output Parameters:
4264: . iscoloring - coloring generated using colorarray information
4266: Level: developer
4268: .seealso: MatGetRowIJ(), MatGetColumnIJ()
4270: @*/
4271: int MatColoringPatch(Mat mat,int n,int ncolors,int *colorarray,ISColoring *iscoloring)
4272: {
4278: MatPreallocated(mat);
4281: if (!mat->ops->coloringpatch){
4282: ISColoringCreate(mat->comm,n,colorarray,iscoloring);
4283: } else {
4284: (*mat->ops->coloringpatch)(mat,n,ncolors,colorarray,iscoloring);
4285: }
4286: return(0);
4287: }
4290: #undef __FUNCT__
4292: /*@
4293: MatSetUnfactored - Resets a factored matrix to be treated as unfactored.
4295: Collective on Mat
4297: Input Parameter:
4298: . mat - the factored matrix to be reset
4300: Notes:
4301: This routine should be used only with factored matrices formed by in-place
4302: factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE
4303: format). This option can save memory, for example, when solving nonlinear
4304: systems with a matrix-free Newton-Krylov method and a matrix-based, in-place
4305: ILU(0) preconditioner.
4307: Note that one can specify in-place ILU(0) factorization by calling
4308: .vb
4309: PCType(pc,PCILU);
4310: PCILUSeUseInPlace(pc);
4311: .ve
4312: or by using the options -pc_type ilu -pc_ilu_in_place
4314: In-place factorization ILU(0) can also be used as a local
4315: solver for the blocks within the block Jacobi or additive Schwarz
4316: methods (runtime option: -sub_pc_ilu_in_place). See the discussion
4317: of these preconditioners in the users manual for details on setting
4318: local solver options.
4320: Most users should employ the simplified SLES interface for linear solvers
4321: instead of working directly with matrix algebra routines such as this.
4322: See, e.g., SLESCreate().
4324: Level: developer
4326: .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace()
4328: Concepts: matrices^unfactored
4330: @*/
4331: int MatSetUnfactored(Mat mat)
4332: {
4338: MatPreallocated(mat);
4339: mat->factor = 0;
4340: if (!mat->ops->setunfactored) return(0);
4341: (*mat->ops->setunfactored)(mat);
4342: return(0);
4343: }
4345: /*MC
4346: MatGetArrayF90 - Accesses a matrix array from Fortran90.
4348: Synopsis:
4349: MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
4351: Not collective
4353: Input Parameter:
4354: . x - matrix
4356: Output Parameters:
4357: + xx_v - the Fortran90 pointer to the array
4358: - ierr - error code
4360: Example of Usage:
4361: .vb
4362: PetscScalar, pointer xx_v(:)
4363: ....
4364: call MatGetArrayF90(x,xx_v,ierr)
4365: a = xx_v(3)
4366: call MatRestoreArrayF90(x,xx_v,ierr)
4367: .ve
4369: Notes:
4370: Not yet supported for all F90 compilers
4372: Level: advanced
4374: .seealso: MatRestoreArrayF90(), MatGetArray(), MatRestoreArray()
4376: Concepts: matrices^accessing array
4378: M*/
4380: /*MC
4381: MatRestoreArrayF90 - Restores a matrix array that has been
4382: accessed with MatGetArrayF90().
4384: Synopsis:
4385: MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
4387: Not collective
4389: Input Parameters:
4390: + x - matrix
4391: - xx_v - the Fortran90 pointer to the array
4393: Output Parameter:
4394: . ierr - error code
4396: Example of Usage:
4397: .vb
4398: PetscScalar, pointer xx_v(:)
4399: ....
4400: call MatGetArrayF90(x,xx_v,ierr)
4401: a = xx_v(3)
4402: call MatRestoreArrayF90(x,xx_v,ierr)
4403: .ve
4404:
4405: Notes:
4406: Not yet supported for all F90 compilers
4408: Level: advanced
4410: .seealso: MatGetArrayF90(), MatGetArray(), MatRestoreArray()
4412: M*/
4415: #undef __FUNCT__
4417: /*@
4418: MatGetSubMatrix - Gets a single submatrix on the same number of processors
4419: as the original matrix.
4421: Collective on Mat
4423: Input Parameters:
4424: + mat - the original matrix
4425: . isrow - rows this processor should obtain
4426: . iscol - columns for all processors you wish to keep
4427: . csize - number of columns "local" to this processor (does nothing for sequential
4428: matrices). This should match the result from VecGetLocalSize(x,...) if you
4429: plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE
4430: - cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
4432: Output Parameter:
4433: . newmat - the new submatrix, of the same type as the old
4435: Level: advanced
4437: Notes: the iscol argument MUST be the same on each processor. You might be
4438: able to create the iscol argument with ISAllGather().
4440: The first time this is called you should use a cll of MAT_INITIAL_MATRIX,
4441: the MatGetSubMatrix() routine will create the newmat for you. Any additional calls
4442: to this routine with a mat of the same nonzero structure will reuse the matrix
4443: generated the first time.
4445: Concepts: matrices^submatrices
4447: .seealso: MatGetSubMatrices(), ISAllGather()
4448: @*/
4449: int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat)
4450: {
4451: int ierr, size;
4452: Mat *local;
4456: MatPreallocated(mat);
4457: MPI_Comm_size(mat->comm,&size);
4459: /* if original matrix is on just one processor then use submatrix generated */
4460: if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) {
4461: MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);
4462: return(0);
4463: } else if (!mat->ops->getsubmatrix && size == 1) {
4464: ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);
4465: *newmat = *local;
4466: ierr = PetscFree(local);
4467: return(0);
4468: }
4470: if (!mat->ops->getsubmatrix) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4471: (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);
4472: return(0);
4473: }
4475: #undef __FUNCT__
4477: /*@C
4478: MatGetPetscMaps - Returns the maps associated with the matrix.
4480: Not Collective
4482: Input Parameter:
4483: . mat - the matrix
4485: Output Parameters:
4486: + rmap - the row (right) map
4487: - cmap - the column (left) map
4489: Level: developer
4491: Concepts: maps^getting from matrix
4493: @*/
4494: int MatGetPetscMaps(Mat mat,PetscMap *rmap,PetscMap *cmap)
4495: {
4501: MatPreallocated(mat);
4502: (*mat->ops->getmaps)(mat,rmap,cmap);
4503: return(0);
4504: }
4506: /*
4507: Version that works for all PETSc matrices
4508: */
4509: #undef __FUNCT__
4511: int MatGetPetscMaps_Petsc(Mat mat,PetscMap *rmap,PetscMap *cmap)
4512: {
4514: if (rmap) *rmap = mat->rmap;
4515: if (cmap) *cmap = mat->cmap;
4516: return(0);
4517: }
4519: #undef __FUNCT__
4521: /*@
4522: MatSetStashInitialSize - sets the sizes of the matrix stash, that is
4523: used during the assembly process to store values that belong to
4524: other processors.
4526: Not Collective
4528: Input Parameters:
4529: + mat - the matrix
4530: . size - the initial size of the stash.
4531: - bsize - the initial size of the block-stash(if used).
4533: Options Database Keys:
4534: + -matstash_initial_size <size> or <size0,size1,...sizep-1>
4535: - -matstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>
4537: Level: intermediate
4539: Notes:
4540: The block-stash is used for values set with VecSetValuesBlocked() while
4541: the stash is used for values set with VecSetValues()
4543: Run with the option -log_info and look for output of the form
4544: MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
4545: to determine the appropriate value, MM, to use for size and
4546: MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
4547: to determine the value, BMM to use for bsize
4549: Concepts: stash^setting matrix size
4550: Concepts: matrices^stash
4552: @*/
4553: int MatSetStashInitialSize(Mat mat,int size, int bsize)
4554: {
4560: MatPreallocated(mat);
4561: MatStashSetInitialSize_Private(&mat->stash,size);
4562: MatStashSetInitialSize_Private(&mat->bstash,bsize);
4563: return(0);
4564: }
4566: #undef __FUNCT__
4568: /*@
4569: MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of
4570: the matrix
4572: Collective on Mat
4574: Input Parameters:
4575: + mat - the matrix
4576: . x,y - the vectors
4577: - w - where the result is stored
4579: Level: intermediate
4581: Notes:
4582: w may be the same vector as y.
4584: This allows one to use either the restriction or interpolation (its transpose)
4585: matrix to do the interpolation
4587: Concepts: interpolation
4589: .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
4591: @*/
4592: int MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w)
4593: {
4594: int M,N,ierr;
4598: MatPreallocated(A);
4599: MatGetSize(A,&M,&N);
4600: if (N > M) {
4601: MatMultTransposeAdd(A,x,y,w);
4602: } else {
4603: MatMultAdd(A,x,y,w);
4604: }
4605: return(0);
4606: }
4608: #undef __FUNCT__
4610: /*@
4611: MatInterpolate - y = A*x or A'*x depending on the shape of
4612: the matrix
4614: Collective on Mat
4616: Input Parameters:
4617: + mat - the matrix
4618: - x,y - the vectors
4620: Level: intermediate
4622: Notes:
4623: This allows one to use either the restriction or interpolation (its transpose)
4624: matrix to do the interpolation
4626: Concepts: matrices^interpolation
4628: .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
4630: @*/
4631: int MatInterpolate(Mat A,Vec x,Vec y)
4632: {
4633: int M,N,ierr;
4637: MatPreallocated(A);
4638: MatGetSize(A,&M,&N);
4639: if (N > M) {
4640: MatMultTranspose(A,x,y);
4641: } else {
4642: MatMult(A,x,y);
4643: }
4644: return(0);
4645: }
4647: #undef __FUNCT__
4649: /*@
4650: MatRestrict - y = A*x or A'*x
4652: Collective on Mat
4654: Input Parameters:
4655: + mat - the matrix
4656: - x,y - the vectors
4658: Level: intermediate
4660: Notes:
4661: This allows one to use either the restriction or interpolation (its transpose)
4662: matrix to do the restriction
4664: Concepts: matrices^restriction
4666: .seealso: MatMultAdd(), MatMultTransposeAdd(), MatInterpolate()
4668: @*/
4669: int MatRestrict(Mat A,Vec x,Vec y)
4670: {
4671: int M,N,ierr;
4675: MatPreallocated(A);
4676: MatGetSize(A,&M,&N);
4677: if (N > M) {
4678: MatMult(A,x,y);
4679: } else {
4680: MatMultTranspose(A,x,y);
4681: }
4682: return(0);
4683: }
4685: #undef __FUNCT__
4687: /*@C
4688: MatNullSpaceAttach - attaches a null space to a matrix.
4689: This null space will be removed from the resulting vector whenever
4690: MatMult() is called
4692: Collective on Mat
4694: Input Parameters:
4695: + mat - the matrix
4696: - nullsp - the null space object
4698: Level: developer
4700: Notes:
4701: Overwrites any previous null space that may have been attached
4703: Concepts: null space^attaching to matrix
4705: .seealso: MatCreate(), MatNullSpaceCreate()
4706: @*/
4707: int MatNullSpaceAttach(Mat mat,MatNullSpace nullsp)
4708: {
4714: MatPreallocated(mat);
4717: if (mat->nullsp) {
4718: MatNullSpaceDestroy(mat->nullsp);
4719: }
4720: mat->nullsp = nullsp;
4721: PetscObjectReference((PetscObject)nullsp);
4722: return(0);
4723: }
4725: #undef __FUNCT__
4727: /*@
4728: MatICCFactor - Performs in-place incomplete Cholesky factorization of matrix.
4730: Collective on Mat
4732: Input Parameters:
4733: + mat - the matrix
4734: . row - row/column permutation
4735: . fill - expected fill factor >= 1.0
4736: - level - level of fill, for ICC(k)
4738: Notes:
4739: Probably really in-place only when level of fill is zero, otherwise allocates
4740: new space to store factored matrix and deletes previous memory.
4742: Most users should employ the simplified SLES interface for linear solvers
4743: instead of working directly with matrix algebra routines such as this.
4744: See, e.g., SLESCreate().
4746: Level: developer
4748: Concepts: matrices^incomplete Cholesky factorization
4749: Concepts: Cholesky factorization
4751: .seealso: MatICCFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
4752: @*/
4753: int MatICCFactor(Mat mat,IS row,PetscReal fill,int level)
4754: {
4760: MatPreallocated(mat);
4761: if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
4762: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4763: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4764: if (!mat->ops->iccfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4765: (*mat->ops->iccfactor)(mat,row,fill,level);
4766: return(0);
4767: }
4769: #undef __FUNCT__
4771: /*@
4772: MatSetValuesAdic - Sets values computed with ADIC automatic differentiation into a matrix.
4774: Not Collective
4776: Input Parameters:
4777: + mat - the matrix
4778: - v - the values compute with ADIC
4780: Level: developer
4782: Notes:
4783: Must call MatSetColoring() before using this routine. Also this matrix must already
4784: have its nonzero pattern determined.
4786: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
4787: MatSetValues(), MatSetColoring(), MatSetValuesAdifor()
4788: @*/
4789: int MatSetValuesAdic(Mat mat,void *v)
4790: {
4797: if (!mat->assembled) {
4798: SETERRQ(1,"Matrix must be already assembled");
4799: }
4800: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
4801: if (!mat->ops->setvaluesadic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4802: (*mat->ops->setvaluesadic)(mat,v);
4803: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
4804: MatView_Private(mat);
4805: return(0);
4806: }
4809: #undef __FUNCT__
4811: /*@
4812: MatSetColoring - Sets a coloring used by calls to MatSetValuesAdic()
4814: Not Collective
4816: Input Parameters:
4817: + mat - the matrix
4818: - coloring - the coloring
4820: Level: developer
4822: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
4823: MatSetValues(), MatSetValuesAdic()
4824: @*/
4825: int MatSetColoring(Mat mat,ISColoring coloring)
4826: {
4833: if (!mat->assembled) {
4834: SETERRQ(1,"Matrix must be already assembled");
4835: }
4836: if (!mat->ops->setcoloring) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4837: (*mat->ops->setcoloring)(mat,coloring);
4838: return(0);
4839: }
4841: #undef __FUNCT__
4843: /*@
4844: MatSetValuesAdifor - Sets values computed with automatic differentiation into a matrix.
4846: Not Collective
4848: Input Parameters:
4849: + mat - the matrix
4850: . nl - leading dimension of v
4851: - v - the values compute with ADIFOR
4853: Level: developer
4855: Notes:
4856: Must call MatSetColoring() before using this routine. Also this matrix must already
4857: have its nonzero pattern determined.
4859: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
4860: MatSetValues(), MatSetColoring()
4861: @*/
4862: int MatSetValuesAdifor(Mat mat,int nl,void *v)
4863: {
4870: if (!mat->assembled) {
4871: SETERRQ(1,"Matrix must be already assembled");
4872: }
4873: PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
4874: if (!mat->ops->setvaluesadifor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4875: (*mat->ops->setvaluesadifor)(mat,nl,v);
4876: PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
4877: return(0);
4878: }
4880: EXTERN int MatMPIAIJDiagonalScaleLocal(Mat,Vec);
4881: EXTERN int MatMPIBAIJDiagonalScaleLocal(Mat,Vec);
4883: #undef __FUNCT__
4885: /*@
4886: MatDiagonalScaleLocal - Scales columns of a matrix given the scaling values including the
4887: ghosted ones.
4889: Not Collective
4891: Input Parameters:
4892: + mat - the matrix
4893: - diag = the diagonal values, including ghost ones
4895: Level: developer
4897: Notes: Works only for MPIAIJ and MPIBAIJ matrices
4898:
4899: .seealso: MatDiagonalScale()
4900: @*/
4901: int MatDiagonalScaleLocal(Mat mat,Vec diag)
4902: {
4903: int ierr;
4904: PetscTruth flag;
4911: if (!mat->assembled) {
4912: SETERRQ(1,"Matrix must be already assembled");
4913: }
4914: PetscLogEventBegin(MAT_Scale,mat,0,0,0);
4915: PetscTypeCompare((PetscObject)mat,MATMPIAIJ,&flag);
4916: if (flag) {
4917: MatMPIAIJDiagonalScaleLocal(mat,diag);
4918: } else {
4919: PetscTypeCompare((PetscObject)mat,MATMPIBAIJ,&flag);
4920: if (flag) {
4921: MatMPIBAIJDiagonalScaleLocal(mat,diag);
4922: } else {
4923: int size;
4924: MPI_Comm_size(mat->comm,&size);
4925: if (size == 1) {
4926: int n,m;
4927: VecGetSize(diag,&n);
4928: MatGetSize(mat,0,&m);
4929: if (m == n) {
4930: MatDiagonalScale(mat,0,diag);
4931: } else {
4932: SETERRQ(1,"Only supprted for sequential matrices when no ghost points/periodic conditions");
4933: }
4934: } else {
4935: SETERRQ(1,"Only supported for MPIAIJ and MPIBAIJ parallel matrices");
4936: }
4937: }
4938: }
4939: PetscLogEventEnd(MAT_Scale,mat,0,0,0);
4940: return(0);
4941: }
4943: #undef __FUNCT__
4945: /*@
4946: MatGetInertia - Gets the inertia from a factored matrix
4948: Collective on Mat
4950: Input Parameter:
4951: . mat - the matrix
4953: Output Parameters:
4954: + nneg - number of negative eigenvalues
4955: . nzero - number of zero eigenvalues
4956: - npos - number of positive eigenvalues
4958: Level: advanced
4960: Notes: Matrix must have been factored by MatCholeskyFactor()
4963: @*/
4964: int MatGetInertia(Mat mat,int *nneg,int *nzero,int *npos)
4965: {
4966: int ierr;
4971: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
4972: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Numeric factor mat is not assembled");
4973: if (!mat->ops->getinertia) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4974: (*mat->ops->getinertia)(mat,nneg,nzero,npos);
4975: return(0);
4976: }
4978: /* ----------------------------------------------------------------*/
4979: #undef __FUNCT__
4981: /*@
4982: MatSolves - Solves A x = b, given a factored matrix, for a collection of vectors
4984: Collective on Mat and Vecs
4986: Input Parameters:
4987: + mat - the factored matrix
4988: - b - the right-hand-side vectors
4990: Output Parameter:
4991: . x - the result vectors
4993: Notes:
4994: The vectors b and x cannot be the same. I.e., one cannot
4995: call MatSolves(A,x,x).
4997: Notes:
4998: Most users should employ the simplified SLES interface for linear solvers
4999: instead of working directly with matrix algebra routines such as this.
5000: See, e.g., SLESCreate().
5002: Level: developer
5004: Concepts: matrices^triangular solves
5006: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd(), MatSolve()
5007: @*/
5008: int MatSolves(Mat mat,Vecs b,Vecs x)
5009: {
5015: MatPreallocated(mat);
5016: if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
5017: if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
5018: if (mat->M == 0 && mat->N == 0) return(0);
5020: if (!mat->ops->solves) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
5021: PetscLogEventBegin(MAT_Solves,mat,0,0,0);
5022: (*mat->ops->solves)(mat,b,x);
5023: PetscLogEventEnd(MAT_Solves,mat,0,0,0);
5024: return(0);
5025: }