xref: /petsc/src/mat/interface/matrix.c (revision 95d5f7c29374efcfd2ca44c2fe93981fbc2b4454)
1 /*$Id: matrix.c,v 1.361 2000/04/09 04:35:53 bsmith Exp bsmith $*/
2 
3 /*
4    This is where the abstract matrix operations are defined
5 */
6 
7 #include "src/mat/matimpl.h"        /*I "mat.h" I*/
8 #include "src/vec/vecimpl.h"
9 
10 #undef __FUNC__
11 #define __FUNC__ /*<a name=""></a>*/"MatGetRow"
12 /*@C
13    MatGetRow - Gets a row of a matrix.  You MUST call MatRestoreRow()
14    for each row that you get to ensure that your application does
15    not bleed memory.
16 
17    Not Collective
18 
19    Input Parameters:
20 +  mat - the matrix
21 -  row - the row to get
22 
23    Output Parameters:
24 +  ncols -  the number of nonzeros in the row
25 .  cols - if nonzero, the column numbers
26 -  vals - if nonzero, the values
27 
28    Notes:
29    This routine is provided for people who need to have direct access
30    to the structure of a matrix.  We hope that we provide enough
31    high-level matrix routines that few users will need it.
32 
33    MatGetRow() always returns 0-based column indices, regardless of
34    whether the internal representation is 0-based (default) or 1-based.
35 
36    For better efficiency, set cols and/or vals to PETSC_NULL if you do
37    not wish to extract these quantities.
38 
39    The user can only examine the values extracted with MatGetRow();
40    the values cannot be altered.  To change the matrix entries, one
41    must use MatSetValues().
42 
43    You can only have one call to MatGetRow() outstanding for a particular
44    matrix at a time, per processor. MatGetRow() can only obtained rows
45    associated with the given processor, it cannot get rows from the
46    other processors; for that we suggest using MatGetSubMatrices(), then
47    MatGetRow() on the submatrix. The row indix passed to MatGetRows()
48    is in the global number of rows.
49 
50    Fortran Notes:
51    The calling sequence from Fortran is
52 .vb
53    MatGetRow(matrix,row,ncols,cols,values,ierr)
54          Mat     matrix (input)
55          integer row    (input)
56          integer ncols  (output)
57          integer cols(maxcols) (output)
58          double precision (or double complex) values(maxcols) output
59 .ve
60    where maxcols >= maximum nonzeros in any row of the matrix.
61 
62    Caution:
63    Do not try to change the contents of the output arrays (cols and vals).
64    In some cases, this may corrupt the matrix.
65 
66    Level: advanced
67 
68 .keywords: matrix, row, get, extract
69 
70 .seealso: MatRestoreRow(), MatSetValues(), MatGetValues(), MatGetSubmatrices(), MatGetDiagonal()
71 @*/
72 int MatGetRow(Mat mat,int row,int *ncols,int **cols,Scalar **vals)
73 {
74   int   ierr;
75 
76   PetscFunctionBegin;
77   PetscValidHeaderSpecific(mat,MAT_COOKIE);
78   PetscValidIntPointer(ncols);
79   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
80   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
81   if (!mat->ops->getrow) SETERRQ(PETSC_ERR_SUP,0,"");
82   PLogEventBegin(MAT_GetRow,mat,0,0,0);
83   ierr = (*mat->ops->getrow)(mat,row,ncols,cols,vals);CHKERRQ(ierr);
84   PLogEventEnd(MAT_GetRow,mat,0,0,0);
85   PetscFunctionReturn(0);
86 }
87 
88 #undef __FUNC__
89 #define __FUNC__ /*<a name=""></a>*/"MatRestoreRow"
90 /*@C
91    MatRestoreRow - Frees any temporary space allocated by MatGetRow().
92 
93    Not Collective
94 
95    Input Parameters:
96 +  mat - the matrix
97 .  row - the row to get
98 .  ncols, cols - the number of nonzeros and their columns
99 -  vals - if nonzero the column values
100 
101    Notes:
102    This routine should be called after you have finished examining the entries.
103 
104    Fortran Notes:
105    The calling sequence from Fortran is
106 .vb
107    MatRestoreRow(matrix,row,ncols,cols,values,ierr)
108       Mat     matrix (input)
109       integer row    (input)
110       integer ncols  (output)
111       integer cols(maxcols) (output)
112       double precision (or double complex) values(maxcols) output
113 .ve
114    Where maxcols >= maximum nonzeros in any row of the matrix.
115 
116    In Fortran MatRestoreRow() MUST be called after MatGetRow()
117    before another call to MatGetRow() can be made.
118 
119    Level: advanced
120 
121 .keywords: matrix, row, restore
122 
123 .seealso:  MatGetRow()
124 @*/
125 int MatRestoreRow(Mat mat,int row,int *ncols,int **cols,Scalar **vals)
126 {
127   int ierr;
128 
129   PetscFunctionBegin;
130   PetscValidHeaderSpecific(mat,MAT_COOKIE);
131   PetscValidIntPointer(ncols);
132   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
133   if (!mat->ops->restorerow) PetscFunctionReturn(0);
134   ierr = (*mat->ops->restorerow)(mat,row,ncols,cols,vals);CHKERRQ(ierr);
135   PetscFunctionReturn(0);
136 }
137 
138 #undef __FUNC__
139 #define __FUNC__ /*<a name=""></a>*/"MatView"
140 /*@C
141    MatView - Visualizes a matrix object.
142 
143    Collective on Mat
144 
145    Input Parameters:
146 +  mat - the matrix
147 -  ptr - visualization context
148 
149   Notes:
150   The available visualization contexts include
151 +    VIEWER_STDOUT_SELF - standard output (default)
152 .    VIEWER_STDOUT_WORLD - synchronized standard
153         output where only the first processor opens
154         the file.  All other processors send their
155         data to the first processor to print.
156 -     VIEWER_DRAW_WORLD - graphical display of nonzero structure
157 
158    The user can open alternative visualization contexts with
159 +    ViewerASCIIOpen() - Outputs matrix to a specified file
160 .    ViewerBinaryOpen() - Outputs matrix in binary to a
161          specified file; corresponding input uses MatLoad()
162 .    ViewerDrawOpen() - Outputs nonzero matrix structure to
163          an X window display
164 -    ViewerSocketOpen() - Outputs matrix to Socket viewer.
165          Currently only the sequential dense and AIJ
166          matrix types support the Socket viewer.
167 
168    The user can call ViewerSetFormat() to specify the output
169    format of ASCII printed objects (when using VIEWER_STDOUT_SELF,
170    VIEWER_STDOUT_WORLD and ViewerASCIIOpen).  Available formats include
171 +    VIEWER_FORMAT_ASCII_DEFAULT - default, prints matrix contents
172 .    VIEWER_FORMAT_ASCII_MATLAB - prints matrix contents in Matlab format
173 .    VIEWER_FORMAT_ASCII_DENSE - prints entire matrix including zeros
174 .    VIEWER_FORMAT_ASCII_COMMON - prints matrix contents, using a sparse
175          format common among all matrix types
176 .    VIEWER_FORMAT_ASCII_IMPL - prints matrix contents, using an implementation-specific
177          format (which is in many cases the same as the default)
178 .    VIEWER_FORMAT_ASCII_INFO - prints basic information about the matrix
179          size and structure (not the matrix entries)
180 -    VIEWER_FORMAT_ASCII_INFO_LONG - prints more detailed information about
181          the matrix structure
182 
183    Level: beginner
184 
185 .keywords: matrix, view, visualize, output, print, write, draw
186 
187 .seealso: ViewerSetFormat(), ViewerASCIIOpen(), ViewerDrawOpen(),
188           ViewerSocketOpen(), ViewerBinaryOpen(), MatLoad()
189 @*/
190 int MatView(Mat mat,Viewer viewer)
191 {
192   int        format,ierr,rows,cols;
193   PetscTruth isascii;
194   char       *cstr;
195 
196   PetscFunctionBegin;
197   PetscValidHeaderSpecific(mat,MAT_COOKIE);
198   if (!viewer) viewer = VIEWER_STDOUT_(mat->comm);
199   PetscValidHeaderSpecific(viewer,VIEWER_COOKIE);
200   PetscCheckSameComm(mat,viewer);
201   if (!mat->assembled) SETERRQ(1,1,"Must call MatAssemblyBegin/End() before viewing matrix");
202 
203   ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr);
204   if (isascii) {
205     ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr);
206     if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) {
207       ierr = ViewerASCIIPrintf(viewer,"Matrix Object:\n");CHKERRQ(ierr);
208       ierr = ViewerASCIIPushTab(viewer);CHKERRQ(ierr);
209       ierr = MatGetType(mat,PETSC_NULL,&cstr);CHKERRQ(ierr);
210       ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr);
211       ierr = ViewerASCIIPrintf(viewer,"type=%s, rows=%d, cols=%d\n",cstr,rows,cols);CHKERRQ(ierr);
212       if (mat->ops->getinfo) {
213         MatInfo info;
214         ierr = MatGetInfo(mat,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr);
215         ierr = ViewerASCIIPrintf(viewer,"total: nonzeros=%d, allocated nonzeros=%d\n",
216                           (int)info.nz_used,(int)info.nz_allocated);CHKERRQ(ierr);
217       }
218     }
219   }
220   if (mat->ops->view) {
221     ierr = ViewerASCIIPushTab(viewer);CHKERRQ(ierr);
222     ierr = (*mat->ops->view)(mat,viewer);CHKERRQ(ierr);
223     ierr = ViewerASCIIPopTab(viewer);CHKERRQ(ierr);
224   } else if (!isascii) {
225     SETERRQ1(1,1,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
226   }
227   if (isascii) {
228     ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr);
229     if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) {
230       ierr = ViewerASCIIPopTab(viewer);CHKERRQ(ierr);
231     }
232   }
233   PetscFunctionReturn(0);
234 }
235 
236 #undef __FUNC__
237 #define __FUNC__ /*<a name=""></a>*/"MatScaleSystem"
238 /*@C
239    MatScaleSystem - Scale a vector solution and right hand side to
240    match the scaling of a scaled matrix.
241 
242    Collective on Mat
243 
244    Input Parameter:
245 +  mat - the matrix
246 .  x - solution vector (or PETSC_NULL)
247 +  b - right hand side vector (or PETSC_NULL)
248 
249 
250    Notes:
251    For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
252    internally scaled, so this does nothing. For MPIROWBS it
253    permutes and diagonally scales.
254 
255    The SLES methods automatically call this routine when required
256    (via PCPreSolve()) so it is rarely used directly.
257 
258    Level: Developer
259 
260 .keywords: matrix, scale
261 
262 .seealso: MatUseScaledForm(), MatUnScaleSystem()
263 @*/
264 int MatScaleSystem(Mat mat,Vec x,Vec b)
265 {
266   int ierr;
267 
268   PetscFunctionBegin;
269   PetscValidHeaderSpecific(mat,MAT_COOKIE);
270   if (x) {PetscValidHeaderSpecific(x,VEC_COOKIE);PetscCheckSameComm(mat,x);}
271   if (b) {PetscValidHeaderSpecific(b,VEC_COOKIE);PetscCheckSameComm(mat,b);}
272 
273   if (mat->ops->scalesystem) {
274     ierr = (*mat->ops->scalesystem)(mat,x,b);CHKERRQ(ierr);
275   }
276   PetscFunctionReturn(0);
277 }
278 
279 #undef __FUNC__
280 #define __FUNC__ /*<a name=""></a>*/"MatUnScaleSystem"
281 /*@C
282    MatUnScaleSystem - Unscales a vector solution and right hand side to
283    match the original scaling of a scaled matrix.
284 
285    Collective on Mat
286 
287    Input Parameter:
288 +  mat - the matrix
289 .  x - solution vector (or PETSC_NULL)
290 +  b - right hand side vector (or PETSC_NULL)
291 
292 
293    Notes:
294    For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
295    internally scaled, so this does nothing. For MPIROWBS it
296    permutes and diagonally scales.
297 
298    The SLES methods automatically call this routine when required
299    (via PCPreSolve()) so it is rarely used directly.
300 
301    Level: Developer
302 
303 .keywords: matrix, scale
304 
305 .seealso: MatUseScaledForm(), MatScaleSystem()
306 @*/
307 int MatUnScaleSystem(Mat mat,Vec x,Vec b)
308 {
309   int ierr;
310 
311   PetscFunctionBegin;
312   PetscValidHeaderSpecific(mat,MAT_COOKIE);
313   if (x) {PetscValidHeaderSpecific(x,VEC_COOKIE);PetscCheckSameComm(mat,x);}
314   if (b) {PetscValidHeaderSpecific(b,VEC_COOKIE);PetscCheckSameComm(mat,b);}
315   if (mat->ops->unscalesystem) {
316     ierr = (*mat->ops->unscalesystem)(mat,x,b);CHKERRQ(ierr);
317   }
318   PetscFunctionReturn(0);
319 }
320 
321 #undef __FUNC__
322 #define __FUNC__ /*<a name=""></a>*/"MatUseScaledForm"
323 /*@C
324    MatUseScaledForm - For matrix storage formats that scale the
325    matrix (for example MPIRowBS matrices are diagonally scaled on
326    assembly) indicates matrix operations (MatMult() etc) are
327    applied using the scaled matrix.
328 
329    Collective on Mat
330 
331    Input Parameter:
332 +  mat - the matrix
333 -  scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for
334             applying the original matrix
335 
336    Notes:
337    For scaled matrix formats, applying the original, unscaled matrix
338    will be slightly more expensive
339 
340    Level: Developer
341 
342 .keywords: matrix, scale
343 
344 .seealso: MatScaleSystem(), MatUnScaleSystem()
345 @*/
346 int MatUseScaledForm(Mat mat,PetscTruth scaled)
347 {
348   int ierr;
349 
350   PetscFunctionBegin;
351   PetscValidHeaderSpecific(mat,MAT_COOKIE);
352   if (mat->ops->usescaledform) {
353     ierr = (*mat->ops->usescaledform)(mat,scaled);CHKERRQ(ierr);
354   }
355   PetscFunctionReturn(0);
356 }
357 
358 #undef __FUNC__
359 #define __FUNC__ /*<a name=""></a>*/"MatDestroy"
360 /*@C
361    MatDestroy - Frees space taken by a matrix.
362 
363    Collective on Mat
364 
365    Input Parameter:
366 .  mat - the matrix
367 
368    Level: beginner
369 
370 .keywords: matrix, destroy
371 @*/
372 int MatDestroy(Mat mat)
373 {
374   int ierr;
375 
376   PetscFunctionBegin;
377   PetscValidHeaderSpecific(mat,MAT_COOKIE);
378   if (--mat->refct > 0) PetscFunctionReturn(0);
379 
380   /* if memory was published with AMS then destroy it */
381   ierr = PetscObjectDepublish(mat);CHKERRQ(ierr);
382 
383   ierr = (*mat->ops->destroy)(mat);CHKERRQ(ierr);
384   PetscFunctionReturn(0);
385 }
386 
387 #undef __FUNC__
388 #define __FUNC__ /*<a name=""></a>*/"MatValid"
389 /*@
390    MatValid - Checks whether a matrix object is valid.
391 
392    Collective on Mat
393 
394    Input Parameter:
395 .  m - the matrix to check
396 
397    Output Parameter:
398    flg - flag indicating matrix status, either
399    PETSC_TRUE if matrix is valid, or PETSC_FALSE otherwise.
400 
401    Level: developer
402 
403 .keywords: matrix, valid
404 @*/
405 int MatValid(Mat m,PetscTruth *flg)
406 {
407   PetscFunctionBegin;
408   PetscValidIntPointer(flg);
409   if (!m)                           *flg = PETSC_FALSE;
410   else if (m->cookie != MAT_COOKIE) *flg = PETSC_FALSE;
411   else                              *flg = PETSC_TRUE;
412   PetscFunctionReturn(0);
413 }
414 
415 #undef __FUNC__
416 #define __FUNC__ /*<a name=""></a>*/"MatSetValues"
417 /*@
418    MatSetValues - Inserts or adds a block of values into a matrix.
419    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
420    MUST be called after all calls to MatSetValues() have been completed.
421 
422    Not Collective
423 
424    Input Parameters:
425 +  mat - the matrix
426 .  v - a logically two-dimensional array of values
427 .  m, idxm - the number of rows and their global indices
428 .  n, idxn - the number of columns and their global indices
429 -  addv - either ADD_VALUES or INSERT_VALUES, where
430    ADD_VALUES adds values to any existing entries, and
431    INSERT_VALUES replaces existing entries with new values
432 
433    Notes:
434    By default the values, v, are row-oriented and unsorted.
435    See MatSetOption() for other options.
436 
437    Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES
438    options cannot be mixed without intervening calls to the assembly
439    routines.
440 
441    MatSetValues() uses 0-based row and column numbers in Fortran
442    as well as in C.
443 
444    Negative indices may be passed in idxm and idxn, these rows and columns are
445    simply ignored. This allows easily inserting element stiffness matrices
446    with homogeneous Dirchlet boundary conditions that you don't want represented
447    in the matrix.
448 
449    Efficiency Alert:
450    The routine MatSetValuesBlocked() may offer much better efficiency
451    for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
452 
453    Level: beginner
454 
455 .keywords: matrix, insert, add, set, values
456 
457 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
458 @*/
459 int MatSetValues(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v,InsertMode addv)
460 {
461   int ierr;
462 
463   PetscFunctionBegin;
464   if (!m || !n) PetscFunctionReturn(0); /* no values to insert */
465   PetscValidHeaderSpecific(mat,MAT_COOKIE);
466   PetscValidIntPointer(idxm);
467   PetscValidIntPointer(idxn);
468   PetscValidScalarPointer(v);
469   if (mat->insertmode == NOT_SET_VALUES) {
470     mat->insertmode = addv;
471   }
472 #if defined(PETSC_USE_BOPT_g)
473   else if (mat->insertmode != addv) {
474     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values");
475   }
476   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
477 #endif
478 
479   if (mat->assembled) {
480     mat->was_assembled = PETSC_TRUE;
481     mat->assembled     = PETSC_FALSE;
482   }
483   PLogEventBegin(MAT_SetValues,mat,0,0,0);
484   if (!mat->ops->setvalues) SETERRQ(PETSC_ERR_SUP,1,"Not supported for this matrix type");
485   ierr = (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr);
486   PLogEventEnd(MAT_SetValues,mat,0,0,0);
487   PetscFunctionReturn(0);
488 }
489 
490 #undef __FUNC__
491 #define __FUNC__ /*<a name=""></a>*/"MatSetValuesBlocked"
492 /*@
493    MatSetValuesBlocked - Inserts or adds a block of values into a matrix.
494 
495    Not Collective
496 
497    Input Parameters:
498 +  mat - the matrix
499 .  v - a logically two-dimensional array of values
500 .  m, idxm - the number of block rows and their global block indices
501 .  n, idxn - the number of block columns and their global block indices
502 -  addv - either ADD_VALUES or INSERT_VALUES, where
503    ADD_VALUES adds values to any existing entries, and
504    INSERT_VALUES replaces existing entries with new values
505 
506    Notes:
507    By default the values, v, are row-oriented and unsorted. So the layout of
508    v is the same as for MatSetValues(). See MatSetOption() for other options.
509 
510    Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
511    options cannot be mixed without intervening calls to the assembly
512    routines.
513 
514    MatSetValuesBlocked() uses 0-based row and column numbers in Fortran
515    as well as in C.
516 
517    Negative indices may be passed in idxm and idxn, these rows and columns are
518    simply ignored. This allows easily inserting element stiffness matrices
519    with homogeneous Dirchlet boundary conditions that you don't want represented
520    in the matrix.
521 
522    Each time an entry is set within a sparse matrix via MatSetValues(),
523    internal searching must be done to determine where to place the the
524    data in the matrix storage space.  By instead inserting blocks of
525    entries via MatSetValuesBlocked(), the overhead of matrix assembly is
526    reduced.
527 
528    Restrictions:
529    MatSetValuesBlocked() is currently supported only for the block AIJ
530    matrix format (MATSEQBAIJ and MATMPIBAIJ, which are created via
531    MatCreateSeqBAIJ() and MatCreateMPIBAIJ()).
532 
533    Level: intermediate
534 
535 .keywords: matrix, insert, add, set, values
536 
537 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
538 @*/
539 int MatSetValuesBlocked(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v,InsertMode addv)
540 {
541   int ierr;
542 
543   PetscFunctionBegin;
544   if (!m || !n) PetscFunctionReturn(0); /* no values to insert */
545   PetscValidHeaderSpecific(mat,MAT_COOKIE);
546   PetscValidIntPointer(idxm);
547   PetscValidIntPointer(idxn);
548   PetscValidScalarPointer(v);
549   if (mat->insertmode == NOT_SET_VALUES) {
550     mat->insertmode = addv;
551   }
552 #if defined(PETSC_USE_BOPT_g)
553   else if (mat->insertmode != addv) {
554     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values");
555   }
556   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
557 #endif
558 
559   if (mat->assembled) {
560     mat->was_assembled = PETSC_TRUE;
561     mat->assembled     = PETSC_FALSE;
562   }
563   PLogEventBegin(MAT_SetValues,mat,0,0,0);
564   if (!mat->ops->setvaluesblocked) SETERRQ(PETSC_ERR_SUP,1,"Not supported for this matrix type");
565   ierr = (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr);
566   PLogEventEnd(MAT_SetValues,mat,0,0,0);
567   PetscFunctionReturn(0);
568 }
569 
570 /*MC
571    MatSetValue - Set a single entry into a matrix.
572 
573    Synopsis:
574    void MatSetValue(Mat m,int row,int col,Scalar value,InsertMode mode);
575 
576    Not collective
577 
578    Input Parameters:
579 +  m - the matrix
580 .  row - the row location of the entry
581 .  col - the column location of the entry
582 .  value - the value to insert
583 -  mode - either INSERT_VALUES or ADD_VALUES
584 
585    Notes:
586    For efficiency one should use MatSetValues() and set several or many
587    values simultaneously if possible.
588 
589    Note that VecSetValue() does NOT return an error code (since this
590    is checked internally).
591 
592    Level: beginner
593 
594 .seealso: MatSetValues()
595 M*/
596 
597 #undef __FUNC__
598 #define __FUNC__ /*<a name=""></a>*/"MatGetValues"
599 /*@
600    MatGetValues - Gets a block of values from a matrix.
601 
602    Not Collective; currently only returns a local block
603 
604    Input Parameters:
605 +  mat - the matrix
606 .  v - a logically two-dimensional array for storing the values
607 .  m, idxm - the number of rows and their global indices
608 -  n, idxn - the number of columns and their global indices
609 
610    Notes:
611    The user must allocate space (m*n Scalars) for the values, v.
612    The values, v, are then returned in a row-oriented format,
613    analogous to that used by default in MatSetValues().
614 
615    MatGetValues() uses 0-based row and column numbers in
616    Fortran as well as in C.
617 
618    MatGetValues() requires that the matrix has been assembled
619    with MatAssemblyBegin()/MatAssemblyEnd().  Thus, calls to
620    MatSetValues() and MatGetValues() CANNOT be made in succession
621    without intermediate matrix assembly.
622 
623    Level: advanced
624 
625 .keywords: matrix, get, values
626 
627 .seealso: MatGetRow(), MatGetSubmatrices(), MatSetValues()
628 @*/
629 int MatGetValues(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v)
630 {
631   int ierr;
632 
633   PetscFunctionBegin;
634   PetscValidHeaderSpecific(mat,MAT_COOKIE);
635   PetscValidIntPointer(idxm);
636   PetscValidIntPointer(idxn);
637   PetscValidScalarPointer(v);
638   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
639   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
640   if (!mat->ops->getvalues) SETERRQ(PETSC_ERR_SUP,0,"");
641 
642   PLogEventBegin(MAT_GetValues,mat,0,0,0);
643   ierr = (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);CHKERRQ(ierr);
644   PLogEventEnd(MAT_GetValues,mat,0,0,0);
645   PetscFunctionReturn(0);
646 }
647 
648 #undef __FUNC__
649 #define __FUNC__ /*<a name=""></a>*/"MatSetLocalToGlobalMapping"
650 /*@
651    MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
652    the routine MatSetValuesLocal() to allow users to insert matrix entries
653    using a local (per-processor) numbering.
654 
655    Not Collective
656 
657    Input Parameters:
658 +  x - the matrix
659 -  mapping - mapping created with ISLocalToGlobalMappingCreate()
660              or ISLocalToGlobalMappingCreateIS()
661 
662    Level: intermediate
663 
664 .keywords: matrix, set, values, local, global, mapping
665 
666 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
667 @*/
668 int MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping)
669 {
670   int ierr;
671   PetscFunctionBegin;
672   PetscValidHeaderSpecific(x,MAT_COOKIE);
673   PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE);
674   if (x->mapping) {
675     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Mapping already set for matrix");
676   }
677 
678   x->mapping = mapping;
679   ierr = PetscObjectReference((PetscObject)mapping);CHKERRQ(ierr);
680   PetscFunctionReturn(0);
681 }
682 
683 #undef __FUNC__
684 #define __FUNC__ /*<a name=""></a>*/"MatSetLocalToGlobalMappingBlock"
685 /*@
686    MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use
687    by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
688    entries using a local (per-processor) numbering.
689 
690    Not Collective
691 
692    Input Parameters:
693 +  x - the matrix
694 -  mapping - mapping created with ISLocalToGlobalMappingCreate() or
695              ISLocalToGlobalMappingCreateIS()
696 
697    Level: intermediate
698 
699 .keywords: matrix, set, values, local ordering
700 
701 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
702            MatSetValuesBlocked(), MatSetValuesLocal()
703 @*/
704 int MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping mapping)
705 {
706   int ierr;
707   PetscFunctionBegin;
708   PetscValidHeaderSpecific(x,MAT_COOKIE);
709   PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE);
710   if (x->bmapping) {
711     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Mapping already set for matrix");
712   }
713 
714   x->bmapping = mapping;
715   ierr = PetscObjectReference((PetscObject)mapping);CHKERRQ(ierr);
716   PetscFunctionReturn(0);
717 }
718 
719 #undef __FUNC__
720 #define __FUNC__ /*<a name=""></a>*/"MatSetValuesLocal"
721 /*@
722    MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
723    using a local ordering of the nodes.
724 
725    Not Collective
726 
727    Input Parameters:
728 +  x - the matrix
729 .  nrow, irow - number of rows and their local indices
730 .  ncol, icol - number of columns and their local indices
731 .  y -  a logically two-dimensional array of values
732 -  addv - either INSERT_VALUES or ADD_VALUES, where
733    ADD_VALUES adds values to any existing entries, and
734    INSERT_VALUES replaces existing entries with new values
735 
736    Notes:
737    Before calling MatSetValuesLocal(), the user must first set the
738    local-to-global mapping by calling MatSetLocalToGlobalMapping().
739 
740    Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES
741    options cannot be mixed without intervening calls to the assembly
742    routines.
743 
744    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
745    MUST be called after all calls to MatSetValuesLocal() have been completed.
746 
747    Level: intermediate
748 
749 .keywords: matrix, set, values, local ordering
750 
751 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping()
752 @*/
753 int MatSetValuesLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,Scalar *y,InsertMode addv)
754 {
755   int ierr,irowm[2048],icolm[2048];
756 
757   PetscFunctionBegin;
758   PetscValidHeaderSpecific(mat,MAT_COOKIE);
759   PetscValidIntPointer(irow);
760   PetscValidIntPointer(icol);
761   PetscValidScalarPointer(y);
762 
763   if (mat->insertmode == NOT_SET_VALUES) {
764     mat->insertmode = addv;
765   }
766 #if defined(PETSC_USE_BOPT_g)
767   else if (mat->insertmode != addv) {
768     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values");
769   }
770   if (!mat->mapping) {
771     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Local to global never set with MatSetLocalToGlobalMapping()");
772   }
773   if (nrow > 2048 || ncol > 2048) {
774     SETERRQ2(PETSC_ERR_SUP,0,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
775   }
776   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
777 #endif
778 
779   if (mat->assembled) {
780     mat->was_assembled = PETSC_TRUE;
781     mat->assembled     = PETSC_FALSE;
782   }
783   PLogEventBegin(MAT_SetValues,mat,0,0,0);
784   ierr = ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm);CHKERRQ(ierr);
785   ierr = ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);CHKERRQ(ierr);
786   ierr = (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr);
787   PLogEventEnd(MAT_SetValues,mat,0,0,0);
788   PetscFunctionReturn(0);
789 }
790 
791 #undef __FUNC__
792 #define __FUNC__ /*<a name=""></a>*/"MatSetValuesBlockedLocal"
793 /*@
794    MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
795    using a local ordering of the nodes a block at a time.
796 
797    Not Collective
798 
799    Input Parameters:
800 +  x - the matrix
801 .  nrow, irow - number of rows and their local indices
802 .  ncol, icol - number of columns and their local indices
803 .  y -  a logically two-dimensional array of values
804 -  addv - either INSERT_VALUES or ADD_VALUES, where
805    ADD_VALUES adds values to any existing entries, and
806    INSERT_VALUES replaces existing entries with new values
807 
808    Notes:
809    Before calling MatSetValuesBlockedLocal(), the user must first set the
810    local-to-global mapping by calling MatSetLocalToGlobalMappingBlock(),
811    where the mapping MUST be set for matrix blocks, not for matrix elements.
812 
813    Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
814    options cannot be mixed without intervening calls to the assembly
815    routines.
816 
817    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
818    MUST be called after all calls to MatSetValuesBlockedLocal() have been completed.
819 
820    Level: intermediate
821 
822 .keywords: matrix, set, values, blocked, local
823 
824 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked()
825 @*/
826 int MatSetValuesBlockedLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,Scalar *y,InsertMode addv)
827 {
828   int ierr,irowm[2048],icolm[2048];
829 
830   PetscFunctionBegin;
831   PetscValidHeaderSpecific(mat,MAT_COOKIE);
832   PetscValidIntPointer(irow);
833   PetscValidIntPointer(icol);
834   PetscValidScalarPointer(y);
835   if (mat->insertmode == NOT_SET_VALUES) {
836     mat->insertmode = addv;
837   }
838 #if defined(PETSC_USE_BOPT_g)
839   else if (mat->insertmode != addv) {
840     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values");
841   }
842   if (!mat->bmapping) {
843     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Local to global never set with MatSetLocalToGlobalMappingBlock()");
844   }
845   if (nrow > 2048 || ncol > 2048) {
846     SETERRQ2(PETSC_ERR_SUP,0,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
847   }
848   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
849 #endif
850 
851   if (mat->assembled) {
852     mat->was_assembled = PETSC_TRUE;
853     mat->assembled     = PETSC_FALSE;
854   }
855   PLogEventBegin(MAT_SetValues,mat,0,0,0);
856   ierr = ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm);CHKERRQ(ierr);
857   ierr = ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm);CHKERRQ(ierr);
858   ierr = (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr);
859   PLogEventEnd(MAT_SetValues,mat,0,0,0);
860   PetscFunctionReturn(0);
861 }
862 
863 /* --------------------------------------------------------*/
864 #undef __FUNC__
865 #define __FUNC__ /*<a name=""></a>*/"MatMult"
866 /*@
867    MatMult - Computes the matrix-vector product, y = Ax.
868 
869    Collective on Mat and Vec
870 
871    Input Parameters:
872 +  mat - the matrix
873 -  x   - the vector to be multilplied
874 
875    Output Parameters:
876 .  y - the result
877 
878    Notes:
879    The vectors x and y cannot be the same.  I.e., one cannot
880    call MatMult(A,y,y).
881 
882    Level: beginner
883 
884 .keywords: matrix, multiply, matrix-vector product
885 
886 .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
887 @*/
888 int MatMult(Mat mat,Vec x,Vec y)
889 {
890   int ierr;
891 
892   PetscFunctionBegin;
893   PetscValidHeaderSpecific(mat,MAT_COOKIE);
894   PetscValidHeaderSpecific(x,VEC_COOKIE);
895   PetscValidHeaderSpecific(y,VEC_COOKIE);
896   PetscCheckSameComm(mat,x);
897   PetscCheckSameComm(mat,y);
898   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
899   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
900   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"x and y must be different vectors");
901   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
902   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
903   if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: local dim %d %d",mat->m,y->n);
904 
905   PLogEventBegin(MAT_Mult,mat,x,y,0);
906   ierr = (*mat->ops->mult)(mat,x,y);CHKERRQ(ierr);
907   PLogEventEnd(MAT_Mult,mat,x,y,0);
908 
909   PetscFunctionReturn(0);
910 }
911 
912 #undef __FUNC__
913 #define __FUNC__ /*<a name=""></a>*/"MatMultTranspose"
914 /*@
915    MatMultTranspose - Computes matrix transpose times a vector.
916 
917    Collective on Mat and Vec
918 
919    Input Parameters:
920 +  mat - the matrix
921 -  x   - the vector to be multilplied
922 
923    Output Parameters:
924 .  y - the result
925 
926    Notes:
927    The vectors x and y cannot be the same.  I.e., one cannot
928    call MatMultTranspose(A,y,y).
929 
930    Level: beginner
931 
932 .keywords: matrix, multiply, matrix-vector product, transpose
933 
934 .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd()
935 @*/
936 int MatMultTranspose(Mat mat,Vec x,Vec y)
937 {
938   int ierr;
939 
940   PetscFunctionBegin;
941   PetscValidHeaderSpecific(mat,MAT_COOKIE);
942   PetscValidHeaderSpecific(x,VEC_COOKIE);
943   PetscValidHeaderSpecific(y,VEC_COOKIE);
944   PetscCheckSameComm(mat,x);
945   PetscCheckSameComm(mat,y);
946   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
947   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
948   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"x and y must be different vectors");
949   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
950   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
951 
952   PLogEventBegin(MAT_MultTranspose,mat,x,y,0);
953   ierr = (*mat->ops->multtranspose)(mat,x,y);CHKERRQ(ierr);
954   PLogEventEnd(MAT_MultTranspose,mat,x,y,0);
955   PetscFunctionReturn(0);
956 }
957 
958 #undef __FUNC__
959 #define __FUNC__ /*<a name=""></a>*/"MatMultAdd"
960 /*@
961     MatMultAdd -  Computes v3 = v2 + A * v1.
962 
963     Collective on Mat and Vec
964 
965     Input Parameters:
966 +   mat - the matrix
967 -   v1, v2 - the vectors
968 
969     Output Parameters:
970 .   v3 - the result
971 
972     Notes:
973     The vectors v1 and v3 cannot be the same.  I.e., one cannot
974     call MatMultAdd(A,v1,v2,v1).
975 
976     Level: beginner
977 
978 .keywords: matrix, multiply, matrix-vector product, add
979 
980 .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd()
981 @*/
982 int MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
983 {
984   int ierr;
985 
986   PetscFunctionBegin;
987   PetscValidHeaderSpecific(mat,MAT_COOKIE);
988   PetscValidHeaderSpecific(v1,VEC_COOKIE);
989   PetscValidHeaderSpecific(v2,VEC_COOKIE);
990   PetscValidHeaderSpecific(v3,VEC_COOKIE);
991   PetscCheckSameComm(mat,v1);
992   PetscCheckSameComm(mat,v2);
993   PetscCheckSameComm(mat,v3);
994   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
995   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
996   if (mat->N != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v1: global dim %d %d",mat->N,v1->N);
997   if (mat->M != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: global dim %d %d",mat->M,v2->N);
998   if (mat->M != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: global dim %d %d",mat->M,v3->N);
999   if (mat->m != v3->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: local dim %d %d",mat->m,v3->n);
1000   if (mat->m != v2->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: local dim %d %d",mat->m,v2->n);
1001   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,0,"v1 and v3 must be different vectors");
1002 
1003   PLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);
1004   ierr = (*mat->ops->multadd)(mat,v1,v2,v3);CHKERRQ(ierr);
1005   PLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);
1006   PetscFunctionReturn(0);
1007 }
1008 
1009 #undef __FUNC__
1010 #define __FUNC__ /*<a name=""></a>*/"MatMultTransposeAdd"
1011 /*@
1012    MatMultTransposeAdd - Computes v3 = v2 + A' * v1.
1013 
1014    Collective on Mat and Vec
1015 
1016    Input Parameters:
1017 +  mat - the matrix
1018 -  v1, v2 - the vectors
1019 
1020    Output Parameters:
1021 .  v3 - the result
1022 
1023    Notes:
1024    The vectors v1 and v3 cannot be the same.  I.e., one cannot
1025    call MatMultTransposeAdd(A,v1,v2,v1).
1026 
1027    Level: beginner
1028 
1029 .keywords: matrix, multiply, matrix-vector product, transpose, add
1030 
1031 .seealso: MatMultTranspose(), MatMultAdd(), MatMult()
1032 @*/
1033 int MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1034 {
1035   int ierr;
1036 
1037   PetscFunctionBegin;
1038   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1039   PetscValidHeaderSpecific(v1,VEC_COOKIE);
1040   PetscValidHeaderSpecific(v2,VEC_COOKIE);
1041   PetscValidHeaderSpecific(v3,VEC_COOKIE);
1042   PetscCheckSameComm(mat,v1);
1043   PetscCheckSameComm(mat,v2);
1044   PetscCheckSameComm(mat,v3);
1045   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1046   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1047   if (!mat->ops->multtransposeadd) SETERRQ(PETSC_ERR_SUP,0,"");
1048   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,0,"v1 and v3 must be different vectors");
1049   if (mat->M != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v1: global dim %d %d",mat->M,v1->N);
1050   if (mat->N != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: global dim %d %d",mat->N,v2->N);
1051   if (mat->N != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: global dim %d %d",mat->N,v3->N);
1052 
1053   PLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);
1054   ierr = (*mat->ops->multtransposeadd)(mat,v1,v2,v3);CHKERRQ(ierr);
1055   PLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);
1056   PetscFunctionReturn(0);
1057 }
1058 /* ------------------------------------------------------------*/
1059 #undef __FUNC__
1060 #define __FUNC__ /*<a name=""></a>*/"MatGetInfo"
1061 /*@C
1062    MatGetInfo - Returns information about matrix storage (number of
1063    nonzeros, memory, etc.).
1064 
1065    Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used
1066    as the flag
1067 
1068    Input Parameters:
1069 .  mat - the matrix
1070 
1071    Output Parameters:
1072 +  flag - flag indicating the type of parameters to be returned
1073    (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
1074    MAT_GLOBAL_SUM - sum over all processors)
1075 -  info - matrix information context
1076 
1077    Notes:
1078    The MatInfo context contains a variety of matrix data, including
1079    number of nonzeros allocated and used, number of mallocs during
1080    matrix assembly, etc.  Additional information for factored matrices
1081    is provided (such as the fill ratio, number of mallocs during
1082    factorization, etc.).  Much of this info is printed to STDOUT
1083    when using the runtime options
1084 $       -log_info -mat_view_info
1085 
1086    Example for C/C++ Users:
1087    See the file ${PETSC_DIR}/include/mat.h for a complete list of
1088    data within the MatInfo context.  For example,
1089 .vb
1090       MatInfo info;
1091       Mat     A;
1092       double  mal, nz_a, nz_u;
1093 
1094       MatGetInfo(A,MAT_LOCAL,&info);
1095       mal  = info.mallocs;
1096       nz_a = info.nz_allocated;
1097 .ve
1098 
1099    Example for Fortran Users:
1100    Fortran users should declare info as a double precision
1101    array of dimension MAT_INFO_SIZE, and then extract the parameters
1102    of interest.  See the file ${PETSC_DIR}/include/finclude/mat.h
1103    a complete list of parameter names.
1104 .vb
1105       double  precision info(MAT_INFO_SIZE)
1106       double  precision mal, nz_a
1107       Mat     A
1108       integer ierr
1109 
1110       call MatGetInfo(A,MAT_LOCAL,info,ierr)
1111       mal = info(MAT_INFO_MALLOCS)
1112       nz_a = info(MAT_INFO_NZ_ALLOCATED)
1113 .ve
1114 
1115     Level: intermediate
1116 
1117 .keywords: matrix, get, info, storage, nonzeros, memory, fill
1118 @*/
1119 int MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
1120 {
1121   int ierr;
1122 
1123   PetscFunctionBegin;
1124   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1125   PetscValidPointer(info);
1126   if (!mat->ops->getinfo) SETERRQ(PETSC_ERR_SUP,0,"");
1127   ierr = (*mat->ops->getinfo)(mat,flag,info);CHKERRQ(ierr);
1128   PetscFunctionReturn(0);
1129 }
1130 
1131 /* ----------------------------------------------------------*/
1132 #undef __FUNC__
1133 #define __FUNC__ /*<a name=""></a>*/"MatILUDTFactor"
1134 /*@C
1135    MatILUDTFactor - Performs a drop tolerance ILU factorization.
1136 
1137    Collective on Mat
1138 
1139    Input Parameters:
1140 +  mat - the matrix
1141 .  info - information about the factorization to be done
1142 .  row - row permutation
1143 -  col - column permutation
1144 
1145    Output Parameters:
1146 .  fact - the factored matrix
1147 
1148    Level: developer
1149 
1150    Notes:
1151    Most users should employ the simplified SLES interface for linear solvers
1152    instead of working directly with matrix algebra routines such as this.
1153    See, e.g., SLESCreate().
1154 
1155    This is currently only supported for the SeqAIJ matrix format using code
1156    from Yousef Saad's SPARSEKIT2 package. That code is copyright by Yousef
1157    Saad with the GNU copyright.
1158 
1159 .keywords: matrix, factor, LU, in-place
1160 
1161 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
1162 @*/
1163 int MatILUDTFactor(Mat mat,MatILUInfo *info,IS row,IS col,Mat *fact)
1164 {
1165   int ierr;
1166 
1167   PetscFunctionBegin;
1168   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1169   PetscValidPointer(fact);
1170   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1171   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1172   if (!mat->ops->iludtfactor) SETERRQ(PETSC_ERR_SUP,0,"");
1173 
1174   PLogEventBegin(MAT_ILUFactor,mat,row,col,0);
1175   ierr = (*mat->ops->iludtfactor)(mat,info,row,col,fact);CHKERRQ(ierr);
1176   PLogEventEnd(MAT_ILUFactor,mat,row,col,0);
1177 
1178   PetscFunctionReturn(0);
1179 }
1180 
1181 #undef __FUNC__
1182 #define __FUNC__ /*<a name=""></a>*/"MatLUFactor"
1183 /*@
1184    MatLUFactor - Performs in-place LU factorization of matrix.
1185 
1186    Collective on Mat
1187 
1188    Input Parameters:
1189 +  mat - the matrix
1190 .  row - row permutation
1191 .  col - column permutation
1192 -  f - expected fill as ratio of original fill.
1193 
1194    Notes:
1195    Most users should employ the simplified SLES interface for linear solvers
1196    instead of working directly with matrix algebra routines such as this.
1197    See, e.g., SLESCreate().
1198 
1199    This changes the state of the matrix to a factored matrix; it cannot be used
1200    for example with MatSetValues() unless one first calls MatSetUnfactored().
1201 
1202    Level: developer
1203 
1204 .keywords: matrix, factor, LU, in-place
1205 
1206 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(),
1207           MatGetOrdering(), MatSetUnfactored()
1208 
1209 @*/
1210 int MatLUFactor(Mat mat,IS row,IS col,PetscReal f)
1211 {
1212   int ierr;
1213 
1214   PetscFunctionBegin;
1215   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1216   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1217   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1218   if (!mat->ops->lufactor) SETERRQ(PETSC_ERR_SUP,0,"");
1219 
1220   PLogEventBegin(MAT_LUFactor,mat,row,col,0);
1221   ierr = (*mat->ops->lufactor)(mat,row,col,f);CHKERRQ(ierr);
1222   PLogEventEnd(MAT_LUFactor,mat,row,col,0);
1223   PetscFunctionReturn(0);
1224 }
1225 
1226 #undef __FUNC__
1227 #define __FUNC__ /*<a name=""></a>*/"MatILUFactor"
1228 /*@
1229    MatILUFactor - Performs in-place ILU factorization of matrix.
1230 
1231    Collective on Mat
1232 
1233    Input Parameters:
1234 +  mat - the matrix
1235 .  row - row permutation
1236 .  col - column permutation
1237 -  info - structure containing
1238 $      levels - number of levels of fill.
1239 $      expected fill - as ratio of original fill.
1240 $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
1241                 missing diagonal entries)
1242 
1243    Notes:
1244    Probably really in-place only when level of fill is zero, otherwise allocates
1245    new space to store factored matrix and deletes previous memory.
1246 
1247    Most users should employ the simplified SLES interface for linear solvers
1248    instead of working directly with matrix algebra routines such as this.
1249    See, e.g., SLESCreate().
1250 
1251    Level: developer
1252 
1253 .keywords: matrix, factor, ILU, in-place
1254 
1255 .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
1256 @*/
1257 int MatILUFactor(Mat mat,IS row,IS col,MatILUInfo *info)
1258 {
1259   int ierr;
1260 
1261   PetscFunctionBegin;
1262   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1263   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square");
1264   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1265   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1266   if (!mat->ops->ilufactor) SETERRQ(PETSC_ERR_SUP,0,"");
1267 
1268   PLogEventBegin(MAT_ILUFactor,mat,row,col,0);
1269   ierr = (*mat->ops->ilufactor)(mat,row,col,info);CHKERRQ(ierr);
1270   PLogEventEnd(MAT_ILUFactor,mat,row,col,0);
1271   PetscFunctionReturn(0);
1272 }
1273 
1274 #undef __FUNC__
1275 #define __FUNC__ /*<a name=""></a>*/"MatLUFactorSymbolic"
1276 /*@
1277    MatLUFactorSymbolic - Performs symbolic LU factorization of matrix.
1278    Call this routine before calling MatLUFactorNumeric().
1279 
1280    Collective on Mat
1281 
1282    Input Parameters:
1283 +  mat - the matrix
1284 .  row, col - row and column permutations
1285 -  f - expected fill as ratio of the original number of nonzeros,
1286        for example 3.0; choosing this parameter well can result in
1287        more efficient use of time and space. Run with the option -log_info
1288        to determine an optimal value to use
1289 
1290    Output Parameter:
1291 .  fact - new matrix that has been symbolically factored
1292 
1293    Notes:
1294    See the users manual for additional information about
1295    choosing the fill factor for better efficiency.
1296 
1297    Most users should employ the simplified SLES interface for linear solvers
1298    instead of working directly with matrix algebra routines such as this.
1299    See, e.g., SLESCreate().
1300 
1301    Level: developer
1302 
1303 .keywords: matrix, factor, LU, symbolic, fill
1304 
1305 .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor()
1306 @*/
1307 int MatLUFactorSymbolic(Mat mat,IS row,IS col,PetscReal f,Mat *fact)
1308 {
1309   int ierr;
1310 
1311   PetscFunctionBegin;
1312   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1313   PetscValidPointer(fact);
1314   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1315   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1316   if (!mat->ops->lufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"");
1317 
1318   PLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);
1319   ierr = (*mat->ops->lufactorsymbolic)(mat,row,col,f,fact);CHKERRQ(ierr);
1320   PLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);
1321   PetscFunctionReturn(0);
1322 }
1323 
1324 #undef __FUNC__
1325 #define __FUNC__ /*<a name=""></a>*/"MatLUFactorNumeric"
1326 /*@
1327    MatLUFactorNumeric - Performs numeric LU factorization of a matrix.
1328    Call this routine after first calling MatLUFactorSymbolic().
1329 
1330    Collective on Mat
1331 
1332    Input Parameters:
1333 +  mat - the matrix
1334 -  row, col - row and  column permutations
1335 
1336    Output Parameters:
1337 .  fact - symbolically factored matrix that must have been generated
1338           by MatLUFactorSymbolic()
1339 
1340    Notes:
1341    See MatLUFactor() for in-place factorization.  See
1342    MatCholeskyFactorNumeric() for the symmetric, positive definite case.
1343 
1344    Most users should employ the simplified SLES interface for linear solvers
1345    instead of working directly with matrix algebra routines such as this.
1346    See, e.g., SLESCreate().
1347 
1348    Level: developer
1349 
1350 .keywords: matrix, factor, LU, numeric
1351 
1352 .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
1353 @*/
1354 int MatLUFactorNumeric(Mat mat,Mat *fact)
1355 {
1356   int        ierr;
1357   PetscTruth flg;
1358 
1359   PetscFunctionBegin;
1360   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1361   PetscValidPointer(fact);
1362   PetscValidHeaderSpecific(*fact,MAT_COOKIE);
1363   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1364   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1365     SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dimensions are different %d should = %d %d should = %d",
1366             mat->M,(*fact)->M,mat->N,(*fact)->N);
1367   }
1368   if (!(*fact)->ops->lufactornumeric) SETERRQ(PETSC_ERR_SUP,0,"");
1369 
1370   PLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0);
1371   ierr = (*(*fact)->ops->lufactornumeric)(mat,fact);CHKERRQ(ierr);
1372   PLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0);
1373   ierr = OptionsHasName(PETSC_NULL,"-mat_view_draw",&flg);CHKERRQ(ierr);
1374   if (flg) {
1375     ierr = OptionsHasName(PETSC_NULL,"-mat_view_contour",&flg);CHKERRQ(ierr);
1376     if (flg) {
1377       ViewerPushFormat(VIEWER_DRAW_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr);
1378     }
1379     ierr = MatView(*fact,VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
1380     ierr = ViewerFlush(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
1381     if (flg) {
1382       ViewerPopFormat(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
1383     }
1384   }
1385   PetscFunctionReturn(0);
1386 }
1387 
1388 #undef __FUNC__
1389 #define __FUNC__ /*<a name=""></a>*/"MatCholeskyFactor"
1390 /*@
1391    MatCholeskyFactor - Performs in-place Cholesky factorization of a
1392    symmetric matrix.
1393 
1394    Collective on Mat
1395 
1396    Input Parameters:
1397 +  mat - the matrix
1398 .  perm - row and column permutations
1399 -  f - expected fill as ratio of original fill
1400 
1401    Notes:
1402    See MatLUFactor() for the nonsymmetric case.  See also
1403    MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().
1404 
1405    Most users should employ the simplified SLES interface for linear solvers
1406    instead of working directly with matrix algebra routines such as this.
1407    See, e.g., SLESCreate().
1408 
1409    Level: developer
1410 
1411 .keywords: matrix, factor, in-place, Cholesky
1412 
1413 .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
1414           MatGetOrdering()
1415 
1416 @*/
1417 int MatCholeskyFactor(Mat mat,IS perm,PetscReal f)
1418 {
1419   int ierr;
1420 
1421   PetscFunctionBegin;
1422   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1423   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Matrix must be square");
1424   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1425   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1426   if (!mat->ops->choleskyfactor) SETERRQ(PETSC_ERR_SUP,0,"");
1427 
1428   PLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);
1429   ierr = (*mat->ops->choleskyfactor)(mat,perm,f);CHKERRQ(ierr);
1430   PLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);
1431   PetscFunctionReturn(0);
1432 }
1433 
1434 #undef __FUNC__
1435 #define __FUNC__ /*<a name=""></a>*/"MatCholeskyFactorSymbolic"
1436 /*@
1437    MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
1438    of a symmetric matrix.
1439 
1440    Collective on Mat
1441 
1442    Input Parameters:
1443 +  mat - the matrix
1444 .  perm - row and column permutations
1445 -  f - expected fill as ratio of original
1446 
1447    Output Parameter:
1448 .  fact - the factored matrix
1449 
1450    Notes:
1451    See MatLUFactorSymbolic() for the nonsymmetric case.  See also
1452    MatCholeskyFactor() and MatCholeskyFactorNumeric().
1453 
1454    Most users should employ the simplified SLES interface for linear solvers
1455    instead of working directly with matrix algebra routines such as this.
1456    See, e.g., SLESCreate().
1457 
1458    Level: developer
1459 
1460 .keywords: matrix, factor, factorization, symbolic, Cholesky
1461 
1462 .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
1463           MatGetOrdering()
1464 
1465 @*/
1466 int MatCholeskyFactorSymbolic(Mat mat,IS perm,PetscReal f,Mat *fact)
1467 {
1468   int ierr;
1469 
1470   PetscFunctionBegin;
1471   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1472   PetscValidPointer(fact);
1473   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Matrix must be square");
1474   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1475   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1476   if (!mat->ops->choleskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"");
1477 
1478   PLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
1479   ierr = (*mat->ops->choleskyfactorsymbolic)(mat,perm,f,fact);CHKERRQ(ierr);
1480   PLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
1481   PetscFunctionReturn(0);
1482 }
1483 
1484 #undef __FUNC__
1485 #define __FUNC__ /*<a name=""></a>*/"MatCholeskyFactorNumeric"
1486 /*@
1487    MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
1488    of a symmetric matrix. Call this routine after first calling
1489    MatCholeskyFactorSymbolic().
1490 
1491    Collective on Mat
1492 
1493    Input Parameter:
1494 .  mat - the initial matrix
1495 
1496    Output Parameter:
1497 .  fact - the factored matrix
1498 
1499    Notes:
1500    Most users should employ the simplified SLES interface for linear solvers
1501    instead of working directly with matrix algebra routines such as this.
1502    See, e.g., SLESCreate().
1503 
1504    Level: developer
1505 
1506 .keywords: matrix, factor, numeric, Cholesky
1507 
1508 .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
1509 @*/
1510 int MatCholeskyFactorNumeric(Mat mat,Mat *fact)
1511 {
1512   int ierr;
1513 
1514   PetscFunctionBegin;
1515   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1516   PetscValidPointer(fact);
1517   if (!mat->ops->choleskyfactornumeric) SETERRQ(PETSC_ERR_SUP,0,"");
1518   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1519   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1520     SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dim %d should = %d %d should = %d",
1521             mat->M,(*fact)->M,mat->N,(*fact)->N);
1522   }
1523 
1524   PLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
1525   ierr = (*mat->ops->choleskyfactornumeric)(mat,fact);CHKERRQ(ierr);
1526   PLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
1527   PetscFunctionReturn(0);
1528 }
1529 
1530 /* ----------------------------------------------------------------*/
1531 #undef __FUNC__
1532 #define __FUNC__ /*<a name=""></a>*/"MatSolve"
1533 /*@
1534    MatSolve - Solves A x = b, given a factored matrix.
1535 
1536    Collective on Mat and Vec
1537 
1538    Input Parameters:
1539 +  mat - the factored matrix
1540 -  b - the right-hand-side vector
1541 
1542    Output Parameter:
1543 .  x - the result vector
1544 
1545    Notes:
1546    The vectors b and x cannot be the same.  I.e., one cannot
1547    call MatSolve(A,x,x).
1548 
1549    Notes:
1550    Most users should employ the simplified SLES interface for linear solvers
1551    instead of working directly with matrix algebra routines such as this.
1552    See, e.g., SLESCreate().
1553 
1554    Level: developer
1555 
1556 .keywords: matrix, linear system, solve, LU, Cholesky, triangular solve
1557 
1558 .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd()
1559 @*/
1560 int MatSolve(Mat mat,Vec b,Vec x)
1561 {
1562   int ierr;
1563 
1564   PetscFunctionBegin;
1565   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1566   PetscValidHeaderSpecific(b,VEC_COOKIE);
1567   PetscValidHeaderSpecific(x,VEC_COOKIE);
1568   PetscCheckSameComm(mat,b);
1569   PetscCheckSameComm(mat,x);
1570   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1571   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1572   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1573   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1574   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1575   if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0);
1576 
1577   if (!mat->ops->solve) SETERRQ(PETSC_ERR_SUP,0,"");
1578   PLogEventBegin(MAT_Solve,mat,b,x,0);
1579   ierr = (*mat->ops->solve)(mat,b,x);CHKERRQ(ierr);
1580   PLogEventEnd(MAT_Solve,mat,b,x,0);
1581   PetscFunctionReturn(0);
1582 }
1583 
1584 #undef __FUNC__
1585 #define __FUNC__ /*<a name=""></a>*/"MatForwardSolve"
1586 /* @
1587    MatForwardSolve - Solves L x = b, given a factored matrix, A = LU.
1588 
1589    Collective on Mat and Vec
1590 
1591    Input Parameters:
1592 +  mat - the factored matrix
1593 -  b - the right-hand-side vector
1594 
1595    Output Parameter:
1596 .  x - the result vector
1597 
1598    Notes:
1599    MatSolve() should be used for most applications, as it performs
1600    a forward solve followed by a backward solve.
1601 
1602    The vectors b and x cannot be the same.  I.e., one cannot
1603    call MatForwardSolve(A,x,x).
1604 
1605    Most users should employ the simplified SLES interface for linear solvers
1606    instead of working directly with matrix algebra routines such as this.
1607    See, e.g., SLESCreate().
1608 
1609    Level: developer
1610 
1611 .keywords: matrix, forward, LU, Cholesky, triangular solve
1612 
1613 .seealso: MatSolve(), MatBackwardSolve()
1614 @ */
1615 int MatForwardSolve(Mat mat,Vec b,Vec x)
1616 {
1617   int ierr;
1618 
1619   PetscFunctionBegin;
1620   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1621   PetscValidHeaderSpecific(b,VEC_COOKIE);
1622   PetscValidHeaderSpecific(x,VEC_COOKIE);
1623   PetscCheckSameComm(mat,b);
1624   PetscCheckSameComm(mat,x);
1625   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1626   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1627   if (!mat->ops->forwardsolve) SETERRQ(PETSC_ERR_SUP,0,"");
1628   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1629   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1630   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1631 
1632   PLogEventBegin(MAT_ForwardSolve,mat,b,x,0);
1633   ierr = (*mat->ops->forwardsolve)(mat,b,x);CHKERRQ(ierr);
1634   PLogEventEnd(MAT_ForwardSolve,mat,b,x,0);
1635   PetscFunctionReturn(0);
1636 }
1637 
1638 #undef __FUNC__
1639 #define __FUNC__ /*<a name=""></a>*/"MatBackwardSolve"
1640 /* @
1641    MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.
1642 
1643    Collective on Mat and Vec
1644 
1645    Input Parameters:
1646 +  mat - the factored matrix
1647 -  b - the right-hand-side vector
1648 
1649    Output Parameter:
1650 .  x - the result vector
1651 
1652    Notes:
1653    MatSolve() should be used for most applications, as it performs
1654    a forward solve followed by a backward solve.
1655 
1656    The vectors b and x cannot be the same.  I.e., one cannot
1657    call MatBackwardSolve(A,x,x).
1658 
1659    Most users should employ the simplified SLES interface for linear solvers
1660    instead of working directly with matrix algebra routines such as this.
1661    See, e.g., SLESCreate().
1662 
1663    Level: developer
1664 
1665 .keywords: matrix, backward, LU, Cholesky, triangular solve
1666 
1667 .seealso: MatSolve(), MatForwardSolve()
1668 @ */
1669 int MatBackwardSolve(Mat mat,Vec b,Vec x)
1670 {
1671   int ierr;
1672 
1673   PetscFunctionBegin;
1674   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1675   PetscValidHeaderSpecific(b,VEC_COOKIE);
1676   PetscValidHeaderSpecific(x,VEC_COOKIE);
1677   PetscCheckSameComm(mat,b);
1678   PetscCheckSameComm(mat,x);
1679   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1680   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1681   if (!mat->ops->backwardsolve) SETERRQ(PETSC_ERR_SUP,0,"");
1682   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1683   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1684   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1685 
1686   PLogEventBegin(MAT_BackwardSolve,mat,b,x,0);
1687   ierr = (*mat->ops->backwardsolve)(mat,b,x);CHKERRQ(ierr);
1688   PLogEventEnd(MAT_BackwardSolve,mat,b,x,0);
1689   PetscFunctionReturn(0);
1690 }
1691 
1692 #undef __FUNC__
1693 #define __FUNC__ /*<a name=""></a>*/"MatSolveAdd"
1694 /*@
1695    MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.
1696 
1697    Collective on Mat and Vec
1698 
1699    Input Parameters:
1700 +  mat - the factored matrix
1701 .  b - the right-hand-side vector
1702 -  y - the vector to be added to
1703 
1704    Output Parameter:
1705 .  x - the result vector
1706 
1707    Notes:
1708    The vectors b and x cannot be the same.  I.e., one cannot
1709    call MatSolveAdd(A,x,y,x).
1710 
1711    Most users should employ the simplified SLES interface for linear solvers
1712    instead of working directly with matrix algebra routines such as this.
1713    See, e.g., SLESCreate().
1714 
1715    Level: developer
1716 
1717 .keywords: matrix, linear system, solve, LU, Cholesky, add
1718 
1719 .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd()
1720 @*/
1721 int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
1722 {
1723   Scalar one = 1.0;
1724   Vec    tmp;
1725   int    ierr;
1726 
1727   PetscFunctionBegin;
1728   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1729   PetscValidHeaderSpecific(y,VEC_COOKIE);
1730   PetscValidHeaderSpecific(b,VEC_COOKIE);
1731   PetscValidHeaderSpecific(x,VEC_COOKIE);
1732   PetscCheckSameComm(mat,b);
1733   PetscCheckSameComm(mat,y);
1734   PetscCheckSameComm(mat,x);
1735   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1736   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1737   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1738   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1739   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
1740   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1741   if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim %d %d",x->n,y->n);
1742 
1743   PLogEventBegin(MAT_SolveAdd,mat,b,x,y);
1744   if (mat->ops->solveadd)  {
1745     ierr = (*mat->ops->solveadd)(mat,b,y,x);CHKERRQ(ierr);
1746   } else {
1747     /* do the solve then the add manually */
1748     if (x != y) {
1749       ierr = MatSolve(mat,b,x);CHKERRQ(ierr);
1750       ierr = VecAXPY(&one,y,x);CHKERRQ(ierr);
1751     } else {
1752       ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr);
1753       PLogObjectParent(mat,tmp);
1754       ierr = VecCopy(x,tmp);CHKERRQ(ierr);
1755       ierr = MatSolve(mat,b,x);CHKERRQ(ierr);
1756       ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr);
1757       ierr = VecDestroy(tmp);CHKERRQ(ierr);
1758     }
1759   }
1760   PLogEventEnd(MAT_SolveAdd,mat,b,x,y);
1761   PetscFunctionReturn(0);
1762 }
1763 
1764 #undef __FUNC__
1765 #define __FUNC__ /*<a name=""></a>*/"MatSolveTranspose"
1766 /*@
1767    MatSolveTranspose - Solves A' x = b, given a factored matrix.
1768 
1769    Collective on Mat and Vec
1770 
1771    Input Parameters:
1772 +  mat - the factored matrix
1773 -  b - the right-hand-side vector
1774 
1775    Output Parameter:
1776 .  x - the result vector
1777 
1778    Notes:
1779    The vectors b and x cannot be the same.  I.e., one cannot
1780    call MatSolveTranspose(A,x,x).
1781 
1782    Most users should employ the simplified SLES interface for linear solvers
1783    instead of working directly with matrix algebra routines such as this.
1784    See, e.g., SLESCreate().
1785 
1786    Level: developer
1787 
1788 .keywords: matrix, linear system, solve, LU, Cholesky, transpose
1789 
1790 .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd()
1791 @*/
1792 int MatSolveTranspose(Mat mat,Vec b,Vec x)
1793 {
1794   int ierr;
1795 
1796   PetscFunctionBegin;
1797   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1798   PetscValidHeaderSpecific(b,VEC_COOKIE);
1799   PetscValidHeaderSpecific(x,VEC_COOKIE);
1800   PetscCheckSameComm(mat,b);
1801   PetscCheckSameComm(mat,x);
1802   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1803   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1804   if (!mat->ops->solvetranspose) SETERRQ(PETSC_ERR_SUP,0,"");
1805   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
1806   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);
1807 
1808   PLogEventBegin(MAT_SolveTranspose,mat,b,x,0);
1809   ierr = (*mat->ops->solvetranspose)(mat,b,x);CHKERRQ(ierr);
1810   PLogEventEnd(MAT_SolveTranspose,mat,b,x,0);
1811   PetscFunctionReturn(0);
1812 }
1813 
1814 #undef __FUNC__
1815 #define __FUNC__ /*<a name=""></a>*/"MatSolveTransposeAdd"
1816 /*@
1817    MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a
1818                       factored matrix.
1819 
1820    Collective on Mat and Vec
1821 
1822    Input Parameters:
1823 +  mat - the factored matrix
1824 .  b - the right-hand-side vector
1825 -  y - the vector to be added to
1826 
1827    Output Parameter:
1828 .  x - the result vector
1829 
1830    Notes:
1831    The vectors b and x cannot be the same.  I.e., one cannot
1832    call MatSolveTransposeAdd(A,x,y,x).
1833 
1834    Most users should employ the simplified SLES interface for linear solvers
1835    instead of working directly with matrix algebra routines such as this.
1836    See, e.g., SLESCreate().
1837 
1838    Level: developer
1839 
1840 .keywords: matrix, linear system, solve, LU, Cholesky, transpose, add
1841 
1842 .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose()
1843 @*/
1844 int MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x)
1845 {
1846   Scalar one = 1.0;
1847   int    ierr;
1848   Vec    tmp;
1849 
1850   PetscFunctionBegin;
1851   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1852   PetscValidHeaderSpecific(y,VEC_COOKIE);
1853   PetscValidHeaderSpecific(b,VEC_COOKIE);
1854   PetscValidHeaderSpecific(x,VEC_COOKIE);
1855   PetscCheckSameComm(mat,b);
1856   PetscCheckSameComm(mat,y);
1857   PetscCheckSameComm(mat,x);
1858   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1859   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1860   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
1861   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);
1862   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
1863   if (x->n != y->n)   SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim %d %d",x->n,y->n);
1864 
1865   PLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);
1866   if (mat->ops->solvetransposeadd) {
1867     ierr = (*mat->ops->solvetransposeadd)(mat,b,y,x);CHKERRQ(ierr);
1868   } else {
1869     /* do the solve then the add manually */
1870     if (x != y) {
1871       ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr);
1872       ierr = VecAXPY(&one,y,x);CHKERRQ(ierr);
1873     } else {
1874       ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr);
1875       PLogObjectParent(mat,tmp);
1876       ierr = VecCopy(x,tmp);CHKERRQ(ierr);
1877       ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr);
1878       ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr);
1879       ierr = VecDestroy(tmp);CHKERRQ(ierr);
1880     }
1881   }
1882   PLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);
1883   PetscFunctionReturn(0);
1884 }
1885 /* ----------------------------------------------------------------*/
1886 
1887 #undef __FUNC__
1888 #define __FUNC__ /*<a name=""></a>*/"MatRelax"
1889 /*@
1890    MatRelax - Computes one relaxation sweep.
1891 
1892    Collective on Mat and Vec
1893 
1894    Input Parameters:
1895 +  mat - the matrix
1896 .  b - the right hand side
1897 .  omega - the relaxation factor
1898 .  flag - flag indicating the type of SOR (see below)
1899 .  shift -  diagonal shift
1900 -  its - the number of iterations
1901 
1902    Output Parameters:
1903 .  x - the solution (can contain an initial guess)
1904 
1905    SOR Flags:
1906 .     SOR_FORWARD_SWEEP - forward SOR
1907 .     SOR_BACKWARD_SWEEP - backward SOR
1908 .     SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
1909 .     SOR_LOCAL_FORWARD_SWEEP - local forward SOR
1910 .     SOR_LOCAL_BACKWARD_SWEEP - local forward SOR
1911 .     SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
1912 .     SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies
1913          upper/lower triangular part of matrix to
1914          vector (with omega)
1915 .     SOR_ZERO_INITIAL_GUESS - zero initial guess
1916 
1917    Notes:
1918    SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
1919    SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings
1920    on each processor.
1921 
1922    Application programmers will not generally use MatRelax() directly,
1923    but instead will employ the SLES/PC interface.
1924 
1925    Notes for Advanced Users:
1926    The flags are implemented as bitwise inclusive or operations.
1927    For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
1928    to specify a zero initial guess for SSOR.
1929 
1930    Most users should employ the simplified SLES interface for linear solvers
1931    instead of working directly with matrix algebra routines such as this.
1932    See, e.g., SLESCreate().
1933 
1934    Level: developer
1935 
1936 .keywords: matrix, relax, relaxation, sweep
1937 @*/
1938 int MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,int its,Vec x)
1939 {
1940   int ierr;
1941 
1942   PetscFunctionBegin;
1943   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1944   PetscValidHeaderSpecific(b,VEC_COOKIE);
1945   PetscValidHeaderSpecific(x,VEC_COOKIE);
1946   PetscCheckSameComm(mat,b);
1947   PetscCheckSameComm(mat,x);
1948   if (!mat->ops->relax) SETERRQ(PETSC_ERR_SUP,0,"");
1949   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1950   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1951   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1952   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1953   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1954 
1955   PLogEventBegin(MAT_Relax,mat,b,x,0);
1956   ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,x);CHKERRQ(ierr);
1957   PLogEventEnd(MAT_Relax,mat,b,x,0);
1958   PetscFunctionReturn(0);
1959 }
1960 
1961 #undef __FUNC__
1962 #define __FUNC__ /*<a name=""></a>*/"MatCopy_Basic"
1963 /*
1964       Default matrix copy routine.
1965 */
1966 int MatCopy_Basic(Mat A,Mat B,MatStructure str)
1967 {
1968   int    ierr,i,rstart,rend,nz,*cwork;
1969   Scalar *vwork;
1970 
1971   PetscFunctionBegin;
1972   ierr = MatZeroEntries(B);CHKERRQ(ierr);
1973   ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
1974   for (i=rstart; i<rend; i++) {
1975     ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr);
1976     ierr = MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);CHKERRQ(ierr);
1977     ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr);
1978   }
1979   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1980   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1981   PetscFunctionReturn(0);
1982 }
1983 
1984 #undef __FUNC__
1985 #define __FUNC__ /*<a name=""></a>*/"MatCopy"
1986 /*@C
1987    MatCopy - Copys a matrix to another matrix.
1988 
1989    Collective on Mat
1990 
1991    Input Parameters:
1992 +  A - the matrix
1993 -  str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN
1994 
1995    Output Parameter:
1996 .  B - where the copy is put
1997 
1998    Notes:
1999    If you use SAME_NONZERO_PATTERN then the zero matrices had better have the
2000    same nonzero pattern or the routine will crash.
2001 
2002    MatCopy() copies the matrix entries of a matrix to another existing
2003    matrix (after first zeroing the second matrix).  A related routine is
2004    MatConvert(), which first creates a new matrix and then copies the data.
2005 
2006    Level: intermediate
2007 
2008 .keywords: matrix, copy, convert
2009 
2010 .seealso: MatConvert()
2011 @*/
2012 int MatCopy(Mat A,Mat B,MatStructure str)
2013 {
2014   int ierr;
2015 
2016   PetscFunctionBegin;
2017   PetscValidHeaderSpecific(A,MAT_COOKIE);
2018   PetscValidHeaderSpecific(B,MAT_COOKIE);
2019   PetscCheckSameComm(A,B);
2020   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2021   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2022   if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat A,Mat B: global dim %d %d",A->M,B->M,
2023                                              A->N,B->N);
2024 
2025   PLogEventBegin(MAT_Copy,A,B,0,0);
2026   if (A->ops->copy) {
2027     ierr = (*A->ops->copy)(A,B,str);CHKERRQ(ierr);
2028   } else { /* generic conversion */
2029     ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr);
2030   }
2031   PLogEventEnd(MAT_Copy,A,B,0,0);
2032   PetscFunctionReturn(0);
2033 }
2034 
2035 static int MatConvertersSet = 0;
2036 static int (*MatConverters[MAX_MATRIX_TYPES][MAX_MATRIX_TYPES])(Mat,MatType,Mat*) =
2037            {{0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
2038             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
2039             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
2040             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
2041             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
2042             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}};
2043 
2044 #undef __FUNC__
2045 #define __FUNC__ /*<a name=""></a>*/"MatConvertRegister"
2046 /*@C
2047     MatConvertRegister - Allows one to register a routine that converts between
2048     two matrix types.
2049 
2050     Not Collective
2051 
2052     Input Parameters:
2053 +   intype - the type of matrix (defined in include/mat.h), for example, MATSEQAIJ.
2054 -   outtype - new matrix type, or MATSAME
2055 
2056     Level: advanced
2057 
2058 .seealso: MatConvertRegisterAll()
2059 @*/
2060 int MatConvertRegister(MatType intype,MatType outtype,int (*converter)(Mat,MatType,Mat*))
2061 {
2062   PetscFunctionBegin;
2063   MatConverters[intype][outtype] = converter;
2064   MatConvertersSet               = 1;
2065   PetscFunctionReturn(0);
2066 }
2067 
2068 #undef __FUNC__
2069 #define __FUNC__ /*<a name=""></a>*/"MatConvert"
2070 /*@C
2071    MatConvert - Converts a matrix to another matrix, either of the same
2072    or different type.
2073 
2074    Collective on Mat
2075 
2076    Input Parameters:
2077 +  mat - the matrix
2078 -  newtype - new matrix type.  Use MATSAME to create a new matrix of the
2079    same type as the original matrix.
2080 
2081    Output Parameter:
2082 .  M - pointer to place new matrix
2083 
2084    Notes:
2085    MatConvert() first creates a new matrix and then copies the data from
2086    the first matrix.  A related routine is MatCopy(), which copies the matrix
2087    entries of one matrix to another already existing matrix context.
2088 
2089    Level: intermediate
2090 
2091 .keywords: matrix, copy, convert
2092 
2093 .seealso: MatCopy(), MatDuplicate()
2094 @*/
2095 int MatConvert(Mat mat,MatType newtype,Mat *M)
2096 {
2097   int ierr;
2098 
2099   PetscFunctionBegin;
2100   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2101   PetscValidPointer(M);
2102   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2103   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2104 
2105   if (newtype > MAX_MATRIX_TYPES || newtype < -1) {
2106     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,1,"Not a valid matrix type");
2107   }
2108   *M  = 0;
2109 
2110   if (!MatConvertersSet) {
2111     ierr = MatLoadRegisterAll();CHKERRQ(ierr);
2112   }
2113 
2114   PLogEventBegin(MAT_Convert,mat,0,0,0);
2115   if ((newtype == mat->type || newtype == MATSAME) && mat->ops->duplicate) {
2116     ierr = (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);CHKERRQ(ierr);
2117   } else {
2118     if (!MatConvertersSet) {
2119       ierr = MatConvertRegisterAll();CHKERRQ(ierr);
2120     }
2121     if (!MatConverters[mat->type][newtype]) {
2122       SETERRQ(PETSC_ERR_ARG_WRONG,1,"Invalid matrix type, or matrix converter not registered");
2123     }
2124     ierr = (*MatConverters[mat->type][newtype])(mat,newtype,M);CHKERRQ(ierr);
2125   }
2126   PLogEventEnd(MAT_Convert,mat,0,0,0);
2127   PetscFunctionReturn(0);
2128 }
2129 
2130 #undef __FUNC__
2131 #define __FUNC__ /*<a name=""></a>*/"MatDuplicate"
2132 /*@C
2133    MatDuplicate - Duplicates a matrix including the non-zero structure.
2134 
2135    Collective on Mat
2136 
2137    Input Parameters:
2138 +  mat - the matrix
2139 -  op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero
2140         values as well or not
2141 
2142    Output Parameter:
2143 .  M - pointer to place new matrix
2144 
2145    Level: intermediate
2146 
2147 .keywords: matrix, copy, convert, duplicate
2148 
2149 .seealso: MatCopy(), MatConvert()
2150 @*/
2151 int MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
2152 {
2153   int ierr;
2154 
2155   PetscFunctionBegin;
2156   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2157   PetscValidPointer(M);
2158   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2159   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2160 
2161   *M  = 0;
2162   PLogEventBegin(MAT_Convert,mat,0,0,0);
2163   if (!mat->ops->duplicate) {
2164     SETERRQ(PETSC_ERR_SUP,1,"Not written for this matrix type");
2165   }
2166   ierr = (*mat->ops->duplicate)(mat,op,M);CHKERRQ(ierr);
2167   PLogEventEnd(MAT_Convert,mat,0,0,0);
2168   PetscFunctionReturn(0);
2169 }
2170 
2171 #undef __FUNC__
2172 #define __FUNC__ /*<a name=""></a>*/"MatGetDiagonal"
2173 /*@
2174    MatGetDiagonal - Gets the diagonal of a matrix.
2175 
2176    Collective on Mat and Vec
2177 
2178    Input Parameters:
2179 +  mat - the matrix
2180 -  v - the vector for storing the diagonal
2181 
2182    Output Parameter:
2183 .  v - the diagonal of the matrix
2184 
2185    Notes:
2186    For the SeqAIJ matrix format, this routine may also be called
2187    on a LU factored matrix; in that case it routines the reciprocal of
2188    the diagonal entries in U. It returns the entries permuted by the
2189    row and column permutation used during the symbolic factorization.
2190 
2191    Level: intermediate
2192 
2193 .keywords: matrix, get, diagonal
2194 
2195 .seealso: MatGetRow(), MatGetSubmatrices(), MatGetSubmatrix()
2196 @*/
2197 int MatGetDiagonal(Mat mat,Vec v)
2198 {
2199   int ierr;
2200 
2201   PetscFunctionBegin;
2202   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2203   PetscValidHeaderSpecific(v,VEC_COOKIE);
2204   /* PetscCheckSameComm(mat,v); Could be MPI vector but Seq matrix cause of two submatrix storage */
2205   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2206   if (!mat->ops->getdiagonal) SETERRQ(PETSC_ERR_SUP,0,"");
2207 
2208   ierr = (*mat->ops->getdiagonal)(mat,v);CHKERRQ(ierr);
2209   PetscFunctionReturn(0);
2210 }
2211 
2212 #undef __FUNC__
2213 #define __FUNC__ /*<a name=""></a>*/"MatTranspose"
2214 /*@C
2215    MatTranspose - Computes an in-place or out-of-place transpose of a matrix.
2216 
2217    Collective on Mat
2218 
2219    Input Parameter:
2220 .  mat - the matrix to transpose
2221 
2222    Output Parameters:
2223 .  B - the transpose (or pass in PETSC_NULL for an in-place transpose)
2224 
2225    Level: intermediate
2226 
2227 .keywords: matrix, transpose
2228 
2229 .seealso: MatMultTranspose(), MatMultTransposeAdd()
2230 @*/
2231 int MatTranspose(Mat mat,Mat *B)
2232 {
2233   int ierr;
2234 
2235   PetscFunctionBegin;
2236   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2237   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2238   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2239   if (!mat->ops->transpose) SETERRQ(PETSC_ERR_SUP,0,"");
2240   ierr = (*mat->ops->transpose)(mat,B);CHKERRQ(ierr);
2241   PetscFunctionReturn(0);
2242 }
2243 
2244 #undef __FUNC__
2245 #define __FUNC__ /*<a name=""></a>*/"MatPermute"
2246 /*@C
2247    MatPermute - Creates a new matrix with rows and columns permuted from the
2248    original.
2249 
2250    Collective on Mat
2251 
2252    Input Parameters:
2253 +  mat - the matrix to permute
2254 .  row - row permutation, each processor supplies only the permutation for its rows
2255 -  col - column permutation, each processor needs the entire column permutation, that is
2256          this is the same size as the total number of columns in the matrix
2257 
2258    Output Parameters:
2259 .  B - the permuted matrix
2260 
2261    Level: advanced
2262 
2263 .keywords: matrix, transpose
2264 
2265 .seealso: MatGetOrdering()
2266 @*/
2267 int MatPermute(Mat mat,IS row,IS col,Mat *B)
2268 {
2269   int ierr;
2270 
2271   PetscFunctionBegin;
2272   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2273   PetscValidHeaderSpecific(row,IS_COOKIE);
2274   PetscValidHeaderSpecific(col,IS_COOKIE);
2275   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2276   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2277   if (!mat->ops->permute) SETERRQ(PETSC_ERR_SUP,0,"");
2278   ierr = (*mat->ops->permute)(mat,row,col,B);CHKERRQ(ierr);
2279   PetscFunctionReturn(0);
2280 }
2281 
2282 #undef __FUNC__
2283 #define __FUNC__ /*<a name=""></a>*/"MatEqual"
2284 /*@
2285    MatEqual - Compares two matrices.
2286 
2287    Collective on Mat
2288 
2289    Input Parameters:
2290 +  A - the first matrix
2291 -  B - the second matrix
2292 
2293    Output Parameter:
2294 .  flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.
2295 
2296    Level: intermediate
2297 
2298 .keywords: matrix, equal, equivalent
2299 @*/
2300 int MatEqual(Mat A,Mat B,PetscTruth *flg)
2301 {
2302   int ierr;
2303 
2304   PetscFunctionBegin;
2305   PetscValidHeaderSpecific(A,MAT_COOKIE);
2306   PetscValidHeaderSpecific(B,MAT_COOKIE);
2307   PetscValidIntPointer(flg);
2308   PetscCheckSameComm(A,B);
2309   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2310   if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2311   if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat A,Mat B: global dim %d %d %d %d",
2312                                              A->M,B->M,A->N,B->N);
2313   if (!A->ops->equal) SETERRQ(PETSC_ERR_SUP,0,"");
2314   ierr = (*A->ops->equal)(A,B,flg);CHKERRQ(ierr);
2315   PetscFunctionReturn(0);
2316 }
2317 
2318 #undef __FUNC__
2319 #define __FUNC__ /*<a name=""></a>*/"MatDiagonalScale"
2320 /*@
2321    MatDiagonalScale - Scales a matrix on the left and right by diagonal
2322    matrices that are stored as vectors.  Either of the two scaling
2323    matrices can be PETSC_NULL.
2324 
2325    Collective on Mat
2326 
2327    Input Parameters:
2328 +  mat - the matrix to be scaled
2329 .  l - the left scaling vector (or PETSC_NULL)
2330 -  r - the right scaling vector (or PETSC_NULL)
2331 
2332    Notes:
2333    MatDiagonalScale() computes A = LAR, where
2334    L = a diagonal matrix, R = a diagonal matrix
2335 
2336    Level: intermediate
2337 
2338 .keywords: matrix, diagonal, scale
2339 
2340 .seealso: MatScale()
2341 @*/
2342 int MatDiagonalScale(Mat mat,Vec l,Vec r)
2343 {
2344   int ierr;
2345 
2346   PetscFunctionBegin;
2347   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2348   if (!mat->ops->diagonalscale) SETERRQ(PETSC_ERR_SUP,0,"");
2349   if (l) {PetscValidHeaderSpecific(l,VEC_COOKIE);PetscCheckSameComm(mat,l);}
2350   if (r) {PetscValidHeaderSpecific(r,VEC_COOKIE);PetscCheckSameComm(mat,r);}
2351   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2352   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2353 
2354   PLogEventBegin(MAT_Scale,mat,0,0,0);
2355   ierr = (*mat->ops->diagonalscale)(mat,l,r);CHKERRQ(ierr);
2356   PLogEventEnd(MAT_Scale,mat,0,0,0);
2357   PetscFunctionReturn(0);
2358 }
2359 
2360 #undef __FUNC__
2361 #define __FUNC__ /*<a name=""></a>*/"MatScale"
2362 /*@
2363     MatScale - Scales all elements of a matrix by a given number.
2364 
2365     Collective on Mat
2366 
2367     Input Parameters:
2368 +   mat - the matrix to be scaled
2369 -   a  - the scaling value
2370 
2371     Output Parameter:
2372 .   mat - the scaled matrix
2373 
2374     Level: intermediate
2375 
2376 .keywords: matrix, scale
2377 
2378 .seealso: MatDiagonalScale()
2379 @*/
2380 int MatScale(Scalar *a,Mat mat)
2381 {
2382   int ierr;
2383 
2384   PetscFunctionBegin;
2385   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2386   PetscValidScalarPointer(a);
2387   if (!mat->ops->scale) SETERRQ(PETSC_ERR_SUP,0,"");
2388   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2389   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2390 
2391   PLogEventBegin(MAT_Scale,mat,0,0,0);
2392   ierr = (*mat->ops->scale)(a,mat);CHKERRQ(ierr);
2393   PLogEventEnd(MAT_Scale,mat,0,0,0);
2394   PetscFunctionReturn(0);
2395 }
2396 
2397 #undef __FUNC__
2398 #define __FUNC__ /*<a name=""></a>*/"MatNorm"
2399 /*@
2400    MatNorm - Calculates various norms of a matrix.
2401 
2402    Collective on Mat
2403 
2404    Input Parameters:
2405 +  mat - the matrix
2406 -  type - the type of norm, NORM_1, NORM_2, NORM_FROBENIUS, NORM_INFINITY
2407 
2408    Output Parameters:
2409 .  norm - the resulting norm
2410 
2411    Level: intermediate
2412 
2413 .keywords: matrix, norm, Frobenius
2414 @*/
2415 int MatNorm(Mat mat,NormType type,PetscReal *norm)
2416 {
2417   int ierr;
2418 
2419   PetscFunctionBegin;
2420   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2421   PetscValidScalarPointer(norm);
2422 
2423   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2424   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2425   if (!mat->ops->norm) SETERRQ(PETSC_ERR_SUP,0,"Not for this matrix type");
2426   ierr = (*mat->ops->norm)(mat,type,norm);CHKERRQ(ierr);
2427   PetscFunctionReturn(0);
2428 }
2429 
2430 /*
2431      This variable is used to prevent counting of MatAssemblyBegin() that
2432    are called from within a MatAssemblyEnd().
2433 */
2434 static int MatAssemblyEnd_InUse = 0;
2435 #undef __FUNC__
2436 #define __FUNC__ /*<a name=""></a>*/"MatAssemblyBegin"
2437 /*@
2438    MatAssemblyBegin - Begins assembling the matrix.  This routine should
2439    be called after completing all calls to MatSetValues().
2440 
2441    Collective on Mat
2442 
2443    Input Parameters:
2444 +  mat - the matrix
2445 -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
2446 
2447    Notes:
2448    MatSetValues() generally caches the values.  The matrix is ready to
2449    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
2450    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
2451    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
2452    using the matrix.
2453 
2454    Level: beginner
2455 
2456 .keywords: matrix, assembly, assemble, begin
2457 
2458 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled()
2459 @*/
2460 int MatAssemblyBegin(Mat mat,MatAssemblyType type)
2461 {
2462   int ierr;
2463 
2464   PetscFunctionBegin;
2465   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2466   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix.\n did you forget to call MatSetUnfactored()?");
2467   if (mat->assembled) {
2468     mat->was_assembled = PETSC_TRUE;
2469     mat->assembled     = PETSC_FALSE;
2470   }
2471   if (!MatAssemblyEnd_InUse) {
2472     PLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);
2473     if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);}
2474     PLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);
2475   } else {
2476     if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);}
2477   }
2478   PetscFunctionReturn(0);
2479 }
2480 
2481 #undef __FUNC__
2482 #define __FUNC__ /*<a name=""></a>*/"MatAssembed"
2483 /*@
2484    MatAssembled - Indicates if a matrix has been assembled and is ready for
2485      use; for example, in matrix-vector product.
2486 
2487    Collective on Mat
2488 
2489    Input Parameter:
2490 .  mat - the matrix
2491 
2492    Output Parameter:
2493 .  assembled - PETSC_TRUE or PETSC_FALSE
2494 
2495    Level: advanced
2496 
2497 .keywords: matrix, assembly, assemble, begin
2498 
2499 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin()
2500 @*/
2501 int MatAssembled(Mat mat,PetscTruth *assembled)
2502 {
2503   PetscFunctionBegin;
2504   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2505   *assembled = mat->assembled;
2506   PetscFunctionReturn(0);
2507 }
2508 
2509 #undef __FUNC__
2510 #define __FUNC__ /*<a name=""></a>*/"MatView_Private"
2511 /*
2512     Processes command line options to determine if/how a matrix
2513   is to be viewed. Called by MatAssemblyEnd() and MatLoad().
2514 */
2515 int MatView_Private(Mat mat)
2516 {
2517   int        ierr;
2518   PetscTruth flg;
2519 
2520   PetscFunctionBegin;
2521   ierr = OptionsHasName(mat->prefix,"-mat_view_info",&flg);CHKERRQ(ierr);
2522   if (flg) {
2523     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO,0);CHKERRQ(ierr);
2524     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2525     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2526   }
2527   ierr = OptionsHasName(mat->prefix,"-mat_view_info_detailed",&flg);CHKERRQ(ierr);
2528   if (flg) {
2529     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO_LONG,0);CHKERRQ(ierr);
2530     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2531     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2532   }
2533   ierr = OptionsHasName(mat->prefix,"-mat_view",&flg);CHKERRQ(ierr);
2534   if (flg) {
2535     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2536   }
2537   ierr = OptionsHasName(mat->prefix,"-mat_view_matlab",&flg);CHKERRQ(ierr);
2538   if (flg) {
2539     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_MATLAB,"M");CHKERRQ(ierr);
2540     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2541     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2542   }
2543   ierr = OptionsHasName(mat->prefix,"-mat_view_draw",&flg);CHKERRQ(ierr);
2544   if (flg) {
2545     ierr = OptionsHasName(mat->prefix,"-mat_view_contour",&flg);CHKERRQ(ierr);
2546     if (flg) {
2547       ViewerPushFormat(VIEWER_DRAW_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr);
2548     }
2549     ierr = MatView(mat,VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
2550     ierr = ViewerFlush(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
2551     if (flg) {
2552       ViewerPopFormat(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
2553     }
2554   }
2555   ierr = OptionsHasName(mat->prefix,"-mat_view_socket",&flg);CHKERRQ(ierr);
2556   if (flg) {
2557     ierr = MatView(mat,VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr);
2558     ierr = ViewerFlush(VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr);
2559   }
2560   PetscFunctionReturn(0);
2561 }
2562 
2563 #undef __FUNC__
2564 #define __FUNC__ /*<a name=""></a>*/"MatAssemblyEnd"
2565 /*@
2566    MatAssemblyEnd - Completes assembling the matrix.  This routine should
2567    be called after MatAssemblyBegin().
2568 
2569    Collective on Mat
2570 
2571    Input Parameters:
2572 +  mat - the matrix
2573 -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
2574 
2575    Options Database Keys:
2576 +  -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
2577 .  -mat_view_info_detailed - Prints more detailed info
2578 .  -mat_view - Prints matrix in ASCII format
2579 .  -mat_view_matlab - Prints matrix in Matlab format
2580 .  -mat_view_draw - Draws nonzero structure of matrix, using MatView() and DrawOpenX().
2581 .  -display <name> - Sets display name (default is host)
2582 -  -draw_pause <sec> - Sets number of seconds to pause after display
2583 
2584    Notes:
2585    MatSetValues() generally caches the values.  The matrix is ready to
2586    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
2587    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
2588    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
2589    using the matrix.
2590 
2591    Level: beginner
2592 
2593 .keywords: matrix, assembly, assemble, end
2594 
2595 .seealso: MatAssemblyBegin(), MatSetValues(), DrawOpenX(), MatView(), MatAssembled()
2596 @*/
2597 int MatAssemblyEnd(Mat mat,MatAssemblyType type)
2598 {
2599   int        ierr;
2600   static int inassm = 0;
2601 
2602   PetscFunctionBegin;
2603   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2604 
2605   inassm++;
2606   MatAssemblyEnd_InUse++;
2607   if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
2608     PLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);
2609     if (mat->ops->assemblyend) {
2610       ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr);
2611     }
2612     PLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);
2613   } else {
2614     if (mat->ops->assemblyend) {
2615       ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr);
2616     }
2617   }
2618 
2619   /* Flush assembly is not a true assembly */
2620   if (type != MAT_FLUSH_ASSEMBLY) {
2621     mat->assembled  = PETSC_TRUE; mat->num_ass++;
2622   }
2623   mat->insertmode = NOT_SET_VALUES;
2624   MatAssemblyEnd_InUse--;
2625 
2626   if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) {
2627     ierr = MatView_Private(mat);CHKERRQ(ierr);
2628   }
2629   inassm--;
2630   PetscFunctionReturn(0);
2631 }
2632 
2633 
2634 #undef __FUNC__
2635 #define __FUNC__ /*<a name=""></a>*/"MatCompress"
2636 /*@
2637    MatCompress - Tries to store the matrix in as little space as
2638    possible.  May fail if memory is already fully used, since it
2639    tries to allocate new space.
2640 
2641    Collective on Mat
2642 
2643    Input Parameters:
2644 .  mat - the matrix
2645 
2646    Level: advanced
2647 
2648 .keywords: matrix, compress
2649 @*/
2650 int MatCompress(Mat mat)
2651 {
2652   int ierr;
2653 
2654   PetscFunctionBegin;
2655   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2656   if (mat->ops->compress) {ierr = (*mat->ops->compress)(mat);CHKERRQ(ierr);}
2657   PetscFunctionReturn(0);
2658 }
2659 
2660 #undef __FUNC__
2661 #define __FUNC__ /*<a name=""></a>*/"MatSetOption"
2662 /*@
2663    MatSetOption - Sets a parameter option for a matrix. Some options
2664    may be specific to certain storage formats.  Some options
2665    determine how values will be inserted (or added). Sorted,
2666    row-oriented input will generally assemble the fastest. The default
2667    is row-oriented, nonsorted input.
2668 
2669    Collective on Mat
2670 
2671    Input Parameters:
2672 +  mat - the matrix
2673 -  option - the option, one of those listed below (and possibly others),
2674              e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR
2675 
2676    Options Describing Matrix Structure:
2677 +    MAT_SYMMETRIC - symmetric in terms of both structure and value
2678 -    MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure
2679 
2680    Options For Use with MatSetValues():
2681    Insert a logically dense subblock, which can be
2682 +    MAT_ROW_ORIENTED - row-oriented
2683 .    MAT_COLUMN_ORIENTED - column-oriented
2684 .    MAT_ROWS_SORTED - sorted by row
2685 .    MAT_ROWS_UNSORTED - not sorted by row
2686 .    MAT_COLUMNS_SORTED - sorted by column
2687 -    MAT_COLUMNS_UNSORTED - not sorted by column
2688 
2689    Not these options reflect the data you pass in with MatSetValues(); it has
2690    nothing to do with how the data is stored internally in the matrix
2691    data structure.
2692 
2693    When (re)assembling a matrix, we can restrict the input for
2694    efficiency/debugging purposes.  These options include
2695 +    MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be
2696         allowed if they generate a new nonzero
2697 .    MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed
2698 .    MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if
2699          they generate a nonzero in a new diagonal (for block diagonal format only)
2700 .    MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only)
2701 .    MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries
2702 .    MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry
2703 -    MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly
2704 
2705    Notes:
2706    Some options are relevant only for particular matrix types and
2707    are thus ignored by others.  Other options are not supported by
2708    certain matrix types and will generate an error message if set.
2709 
2710    If using a Fortran 77 module to compute a matrix, one may need to
2711    use the column-oriented option (or convert to the row-oriented
2712    format).
2713 
2714    MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion
2715    that would generate a new entry in the nonzero structure is instead
2716    ignored.  Thus, if memory has not alredy been allocated for this particular
2717    data, then the insertion is ignored. For dense matrices, in which
2718    the entire array is allocated, no entries are ever ignored.
2719 
2720    MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion
2721    that would generate a new entry in the nonzero structure instead produces
2722    an error. (Currently supported for AIJ and BAIJ formats only.)
2723    This is a useful flag when using SAME_NONZERO_PATTERN in calling
2724    SLESSetOperators() to ensure that the nonzero pattern truely does
2725    remain unchanged.
2726 
2727    MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion
2728    that would generate a new entry that has not been preallocated will
2729    instead produce an error. (Currently supported for AIJ and BAIJ formats
2730    only.) This is a useful flag when debugging matrix memory preallocation.
2731 
2732    MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for
2733    other processors should be dropped, rather than stashed.
2734    This is useful if you know that the "owning" processor is also
2735    always generating the correct matrix entries, so that PETSc need
2736    not transfer duplicate entries generated on another processor.
2737 
2738    MAT_USE_HASH_TABLE indicates that a hash table be used to improve the
2739    searches during matrix assembly. When this flag is set, the hash table
2740    is created during the first Matrix Assembly. This hash table is
2741    used the next time through, during MatSetVaules()/MatSetVaulesBlocked()
2742    to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag
2743    should be used with MAT_USE_HASH_TABLE flag. This option is currently
2744    supported by MATMPIBAIJ format only.
2745 
2746    MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries
2747    are kept in the nonzero structure
2748 
2749    MAT_IGNORE_ZERO_ENTRIES - when using ADD_VALUES for AIJ matrices this will stop
2750    zero values from creating a zero location in the matrix
2751 
2752    Level: intermediate
2753 
2754 .keywords: matrix, option, row-oriented, column-oriented, sorted, nonzero
2755 @*/
2756 int MatSetOption(Mat mat,MatOption op)
2757 {
2758   int ierr;
2759 
2760   PetscFunctionBegin;
2761   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2762   if (mat->ops->setoption) {ierr = (*mat->ops->setoption)(mat,op);CHKERRQ(ierr);}
2763   PetscFunctionReturn(0);
2764 }
2765 
2766 #undef __FUNC__
2767 #define __FUNC__ /*<a name=""></a>*/"MatZeroEntries"
2768 /*@
2769    MatZeroEntries - Zeros all entries of a matrix.  For sparse matrices
2770    this routine retains the old nonzero structure.
2771 
2772    Collective on Mat
2773 
2774    Input Parameters:
2775 .  mat - the matrix
2776 
2777    Level: intermediate
2778 
2779 .keywords: matrix, zero, entries
2780 
2781 .seealso: MatZeroRows()
2782 @*/
2783 int MatZeroEntries(Mat mat)
2784 {
2785   int ierr;
2786 
2787   PetscFunctionBegin;
2788   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2789   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2790   if (!mat->ops->zeroentries) SETERRQ(PETSC_ERR_SUP,0,"");
2791 
2792   PLogEventBegin(MAT_ZeroEntries,mat,0,0,0);
2793   ierr = (*mat->ops->zeroentries)(mat);CHKERRQ(ierr);
2794   PLogEventEnd(MAT_ZeroEntries,mat,0,0,0);
2795   PetscFunctionReturn(0);
2796 }
2797 
2798 #undef __FUNC__
2799 #define __FUNC__ /*<a name=""></a>*/"MatZeroRows"
2800 /*@C
2801    MatZeroRows - Zeros all entries (except possibly the main diagonal)
2802    of a set of rows of a matrix.
2803 
2804    Collective on Mat
2805 
2806    Input Parameters:
2807 +  mat - the matrix
2808 .  is - index set of rows to remove
2809 -  diag - pointer to value put in all diagonals of eliminated rows.
2810           Note that diag is not a pointer to an array, but merely a
2811           pointer to a single value.
2812 
2813    Notes:
2814    For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
2815    but does not release memory.  For the dense and block diagonal
2816    formats this does not alter the nonzero structure.
2817 
2818    If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
2819    of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
2820    merely zeroed.
2821 
2822    The user can set a value in the diagonal entry (or for the AIJ and
2823    row formats can optionally remove the main diagonal entry from the
2824    nonzero structure as well, by passing a null pointer (PETSC_NULL
2825    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
2826 
2827    For the parallel case, all processes that share the matrix (i.e.,
2828    those in the communicator used for matrix creation) MUST call this
2829    routine, regardless of whether any rows being zeroed are owned by
2830    them.
2831 
2832 
2833    Level: intermediate
2834 
2835 .keywords: matrix, zero, rows, boundary conditions
2836 
2837 .seealso: MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
2838 @*/
2839 int MatZeroRows(Mat mat,IS is,Scalar *diag)
2840 {
2841   int ierr;
2842 
2843   PetscFunctionBegin;
2844   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2845   PetscValidHeaderSpecific(is,IS_COOKIE);
2846   if (diag) PetscValidScalarPointer(diag);
2847   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2848   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2849   if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,0,"");
2850 
2851   ierr = (*mat->ops->zerorows)(mat,is,diag);CHKERRQ(ierr);
2852   ierr = MatView_Private(mat);CHKERRQ(ierr);
2853   PetscFunctionReturn(0);
2854 }
2855 
2856 #undef __FUNC__
2857 #define __FUNC__ /*<a name=""></a>*/"MatZeroRowsLocal"
2858 /*@C
2859    MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
2860    of a set of rows of a matrix; using local numbering of rows.
2861 
2862    Collective on Mat
2863 
2864    Input Parameters:
2865 +  mat - the matrix
2866 .  is - index set of rows to remove
2867 -  diag - pointer to value put in all diagonals of eliminated rows.
2868           Note that diag is not a pointer to an array, but merely a
2869           pointer to a single value.
2870 
2871    Notes:
2872    For the AIJ matrix formats this removes the old nonzero structure,
2873    but does not release memory.  For the dense and block diagonal
2874    formats this does not alter the nonzero structure.
2875 
2876    The user can set a value in the diagonal entry (or for the AIJ and
2877    row formats can optionally remove the main diagonal entry from the
2878    nonzero structure as well, by passing a null pointer (PETSC_NULL
2879    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
2880 
2881    Level: intermediate
2882 
2883 .keywords: matrix, zero, rows, boundary conditions
2884 
2885 .seealso: MatZeroEntries(), MatZeroRows()
2886 @*/
2887 int MatZeroRowsLocal(Mat mat,IS is,Scalar *diag)
2888 {
2889   int ierr;
2890   IS  newis;
2891 
2892   PetscFunctionBegin;
2893   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2894   PetscValidHeaderSpecific(is,IS_COOKIE);
2895   if (diag) PetscValidScalarPointer(diag);
2896   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2897   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2898   if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,0,"");
2899 
2900   ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);CHKERRQ(ierr);
2901   ierr = (*mat->ops->zerorows)(mat,newis,diag);CHKERRQ(ierr);
2902   ierr = ISDestroy(newis);CHKERRQ(ierr);
2903   PetscFunctionReturn(0);
2904 }
2905 
2906 #undef __FUNC__
2907 #define __FUNC__ /*<a name=""></a>*/"MatGetSize"
2908 /*@
2909    MatGetSize - Returns the numbers of rows and columns in a matrix.
2910 
2911    Not Collective
2912 
2913    Input Parameter:
2914 .  mat - the matrix
2915 
2916    Output Parameters:
2917 +  m - the number of global rows
2918 -  n - the number of global columns
2919 
2920    Level: beginner
2921 
2922 .keywords: matrix, dimension, size, rows, columns, global, get
2923 
2924 .seealso: MatGetLocalSize()
2925 @*/
2926 int MatGetSize(Mat mat,int *m,int* n)
2927 {
2928   int ierr;
2929 
2930   PetscFunctionBegin;
2931   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2932   ierr = (*mat->ops->getsize)(mat,m,n);CHKERRQ(ierr);
2933   PetscFunctionReturn(0);
2934 }
2935 
2936 #undef __FUNC__
2937 #define __FUNC__ /*<a name=""></a>*/"MatGetLocalSize"
2938 /*@
2939    MatGetLocalSize - Returns the number of rows and columns in a matrix
2940    stored locally.  This information may be implementation dependent, so
2941    use with care.
2942 
2943    Not Collective
2944 
2945    Input Parameters:
2946 .  mat - the matrix
2947 
2948    Output Parameters:
2949 +  m - the number of local rows
2950 -  n - the number of local columns
2951 
2952    Level: beginner
2953 
2954 .keywords: matrix, dimension, size, local, rows, columns, get
2955 
2956 .seealso: MatGetSize()
2957 @*/
2958 int MatGetLocalSize(Mat mat,int *m,int* n)
2959 {
2960   int ierr;
2961 
2962   PetscFunctionBegin;
2963   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2964   ierr = (*mat->ops->getlocalsize)(mat,m,n);CHKERRQ(ierr);
2965   PetscFunctionReturn(0);
2966 }
2967 
2968 #undef __FUNC__
2969 #define __FUNC__ /*<a name=""></a>*/"MatGetOwnershipRange"
2970 /*@
2971    MatGetOwnershipRange - Returns the range of matrix rows owned by
2972    this processor, assuming that the matrix is laid out with the first
2973    n1 rows on the first processor, the next n2 rows on the second, etc.
2974    For certain parallel layouts this range may not be well defined.
2975 
2976    Not Collective
2977 
2978    Input Parameters:
2979 .  mat - the matrix
2980 
2981    Output Parameters:
2982 +  m - the global index of the first local row
2983 -  n - one more than the global index of the last local row
2984 
2985    Level: beginner
2986 
2987 .keywords: matrix, get, range, ownership
2988 @*/
2989 int MatGetOwnershipRange(Mat mat,int *m,int* n)
2990 {
2991   int ierr;
2992 
2993   PetscFunctionBegin;
2994   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2995   if (m) PetscValidIntPointer(m);
2996   if (n) PetscValidIntPointer(n);
2997   if (!mat->ops->getownershiprange) SETERRQ(PETSC_ERR_SUP,0,"");
2998   ierr = (*mat->ops->getownershiprange)(mat,m,n);CHKERRQ(ierr);
2999   PetscFunctionReturn(0);
3000 }
3001 
3002 #undef __FUNC__
3003 #define __FUNC__ /*<a name=""></a>*/"MatILUFactorSymbolic"
3004 /*@
3005    MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
3006    Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric()
3007    to complete the factorization.
3008 
3009    Collective on Mat
3010 
3011    Input Parameters:
3012 +  mat - the matrix
3013 .  row - row permutation
3014 .  column - column permutation
3015 -  info - structure containing
3016 $      levels - number of levels of fill.
3017 $      expected fill - as ratio of original fill.
3018 $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
3019                 missing diagonal entries)
3020 
3021    Output Parameters:
3022 .  fact - new matrix that has been symbolically factored
3023 
3024    Notes:
3025    See the users manual for additional information about
3026    choosing the fill factor for better efficiency.
3027 
3028    Most users should employ the simplified SLES interface for linear solvers
3029    instead of working directly with matrix algebra routines such as this.
3030    See, e.g., SLESCreate().
3031 
3032    Level: developer
3033 
3034 .keywords: matrix, factor, incomplete, ILU, symbolic, fill
3035 
3036 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
3037           MatGetOrdering()
3038 
3039 @*/
3040 int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatILUInfo *info,Mat *fact)
3041 {
3042   int ierr;
3043 
3044   PetscFunctionBegin;
3045   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3046   PetscValidPointer(fact);
3047   if (info && info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Levels of fill negative %d",info->levels);
3048   if (info && info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",info->fill);
3049   if (!mat->ops->ilufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Only MatCreateMPIRowbs() matrices support parallel ILU");
3050   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3051   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3052 
3053   PLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);
3054   ierr = (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr);
3055   PLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);
3056   PetscFunctionReturn(0);
3057 }
3058 
3059 #undef __FUNC__
3060 #define __FUNC__ /*<a name=""></a>*/"MatIncompleteCholeskyFactorSymbolic"
3061 /*@
3062    MatIncompleteCholeskyFactorSymbolic - Performs symbolic incomplete
3063    Cholesky factorization for a symmetric matrix.  Use
3064    MatCholeskyFactorNumeric() to complete the factorization.
3065 
3066    Collective on Mat
3067 
3068    Input Parameters:
3069 +  mat - the matrix
3070 .  perm - row and column permutation
3071 .  fill - levels of fill
3072 -  f - expected fill as ratio of original fill
3073 
3074    Output Parameter:
3075 .  fact - the factored matrix
3076 
3077    Notes:
3078    Currently only no-fill factorization is supported.
3079 
3080    Most users should employ the simplified SLES interface for linear solvers
3081    instead of working directly with matrix algebra routines such as this.
3082    See, e.g., SLESCreate().
3083 
3084    Level: developer
3085 
3086 .keywords: matrix, factor, incomplete, ICC, Cholesky, symbolic, fill
3087 
3088 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor()
3089 @*/
3090 int MatIncompleteCholeskyFactorSymbolic(Mat mat,IS perm,PetscReal f,int fill,Mat *fact)
3091 {
3092   int ierr;
3093 
3094   PetscFunctionBegin;
3095   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3096   PetscValidPointer(fact);
3097   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3098   if (fill < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Fill negative %d",fill);
3099   if (f < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",f);
3100   if (!mat->ops->incompletecholeskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support ICC in parallel");
3101   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3102 
3103   PLogEventBegin(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);
3104   ierr = (*mat->ops->incompletecholeskyfactorsymbolic)(mat,perm,f,fill,fact);CHKERRQ(ierr);
3105   PLogEventEnd(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);
3106   PetscFunctionReturn(0);
3107 }
3108 
3109 #undef __FUNC__
3110 #define __FUNC__ /*<a name=""></a>*/"MatGetArray"
3111 /*@C
3112    MatGetArray - Returns a pointer to the element values in the matrix.
3113    The result of this routine is dependent on the underlying matrix data
3114    structure, and may not even work for certain matrix types.  You MUST
3115    call MatRestoreArray() when you no longer need to access the array.
3116 
3117    Not Collective
3118 
3119    Input Parameter:
3120 .  mat - the matrix
3121 
3122    Output Parameter:
3123 .  v - the location of the values
3124 
3125    Currently returns an array only for the dense formats, giving access to
3126    the local portion of the matrix in the usual Fortran column-oriented format.
3127 
3128    Fortran Note:
3129    This routine is used differently from Fortran, e.g.,
3130 .vb
3131         Mat         mat
3132         Scalar      mat_array(1)
3133         PetscOffset i_mat
3134         int         ierr
3135         call MatGetArray(mat,mat_array,i_mat,ierr)
3136 
3137   C  Access first local entry in matrix; note that array is
3138   C  treated as one dimensional
3139         value = mat_array(i_mat + 1)
3140 
3141         [... other code ...]
3142         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3143 .ve
3144 
3145    See the Fortran chapter of the users manual and
3146    petsc/src/mat/examples/tests for details.
3147 
3148    Level: advanced
3149 
3150 .keywords: matrix, array, elements, values
3151 
3152 .seealso: MatRestoreArray(), MatGetArrayF90()
3153 @*/
3154 int MatGetArray(Mat mat,Scalar **v)
3155 {
3156   int ierr;
3157 
3158   PetscFunctionBegin;
3159   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3160   PetscValidPointer(v);
3161   if (!mat->ops->getarray) SETERRQ(PETSC_ERR_SUP,0,"");
3162   ierr = (*mat->ops->getarray)(mat,v);CHKERRQ(ierr);
3163   PetscFunctionReturn(0);
3164 }
3165 
3166 #undef __FUNC__
3167 #define __FUNC__ /*<a name=""></a>*/"MatRestoreArray"
3168 /*@C
3169    MatRestoreArray - Restores the matrix after MatGetArray() has been called.
3170 
3171    Not Collective
3172 
3173    Input Parameter:
3174 +  mat - the matrix
3175 -  v - the location of the values
3176 
3177    Fortran Note:
3178    This routine is used differently from Fortran, e.g.,
3179 .vb
3180         Mat         mat
3181         Scalar      mat_array(1)
3182         PetscOffset i_mat
3183         int         ierr
3184         call MatGetArray(mat,mat_array,i_mat,ierr)
3185 
3186   C  Access first local entry in matrix; note that array is
3187   C  treated as one dimensional
3188         value = mat_array(i_mat + 1)
3189 
3190         [... other code ...]
3191         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3192 .ve
3193 
3194    See the Fortran chapter of the users manual and
3195    petsc/src/mat/examples/tests for details
3196 
3197    Level: advanced
3198 
3199 .keywords: matrix, array, elements, values, restore
3200 
3201 .seealso: MatGetArray(), MatRestoreArrayF90()
3202 @*/
3203 int MatRestoreArray(Mat mat,Scalar **v)
3204 {
3205   int ierr;
3206 
3207   PetscFunctionBegin;
3208   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3209   PetscValidPointer(v);
3210   if (!mat->ops->restorearray) SETERRQ(PETSC_ERR_SUP,0,"");
3211   ierr = (*mat->ops->restorearray)(mat,v);CHKERRQ(ierr);
3212   PetscFunctionReturn(0);
3213 }
3214 
3215 #undef __FUNC__
3216 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrices"
3217 /*@C
3218    MatGetSubMatrices - Extracts several submatrices from a matrix. If submat
3219    points to an array of valid matrices, they may be reused to store the new
3220    submatrices.
3221 
3222    Collective on Mat
3223 
3224    Input Parameters:
3225 +  mat - the matrix
3226 .  n   - the number of submatrixes to be extracted (on this processor, may be zero)
3227 .  irow, icol - index sets of rows and columns to extract
3228 -  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3229 
3230    Output Parameter:
3231 .  submat - the array of submatrices
3232 
3233    Notes:
3234    MatGetSubMatrices() can extract only sequential submatrices
3235    (from both sequential and parallel matrices). Use MatGetSubMatrix()
3236    to extract a parallel submatrix.
3237 
3238    When extracting submatrices from a parallel matrix, each processor can
3239    form a different submatrix by setting the rows and columns of its
3240    individual index sets according to the local submatrix desired.
3241 
3242    When finished using the submatrices, the user should destroy
3243    them with MatDestroySubMatrices().
3244 
3245    MAT_REUSE_MATRIX can only be used when the nonzero structure of the
3246    original matrix has not changed from that last call to MatGetSubMatrices().
3247 
3248    Fortran Note:
3249    The Fortran interface is slightly different from that given below; it
3250    requires one to pass in  as submat a Mat (integer) array of size at least m.
3251 
3252    Level: advanced
3253 
3254 .keywords: matrix, get, submatrix, submatrices
3255 
3256 .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal()
3257 @*/
3258 int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatReuse scall,Mat **submat)
3259 {
3260   int ierr;
3261 
3262   PetscFunctionBegin;
3263   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3264   if (!mat->ops->getsubmatrices) SETERRQ(PETSC_ERR_SUP,0,"");
3265   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3266 
3267   PLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);
3268   ierr = (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);CHKERRQ(ierr);
3269   PLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);
3270 
3271   PetscFunctionReturn(0);
3272 }
3273 
3274 #undef __FUNC__
3275 #define __FUNC__ /*<a name=""></a>*/"MatDestroyMatrices"
3276 /*@C
3277    MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices().
3278 
3279    Collective on Mat
3280 
3281    Input Parameters:
3282 +  n - the number of local matrices
3283 -  mat - the matrices
3284 
3285    Level: advanced
3286 
3287 .keywords: matrix, destroy, submatrix, submatrices
3288 
3289 .seealso: MatGetSubMatrices()
3290 @*/
3291 int MatDestroyMatrices(int n,Mat **mat)
3292 {
3293   int ierr,i;
3294 
3295   PetscFunctionBegin;
3296   if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,1,"Trying to destroy negative number of matrices %d",n);
3297   PetscValidPointer(mat);
3298   for (i=0; i<n; i++) {
3299     ierr = MatDestroy((*mat)[i]);CHKERRQ(ierr);
3300   }
3301   if (n) {ierr = PetscFree(*mat);CHKERRQ(ierr);}
3302   PetscFunctionReturn(0);
3303 }
3304 
3305 #undef __FUNC__
3306 #define __FUNC__ /*<a name=""></a>*/"MatIncreaseOverlap"
3307 /*@
3308    MatIncreaseOverlap - Given a set of submatrices indicated by index sets,
3309    replaces the index sets by larger ones that represent submatrices with
3310    additional overlap.
3311 
3312    Collective on Mat
3313 
3314    Input Parameters:
3315 +  mat - the matrix
3316 .  n   - the number of index sets
3317 .  is  - the array of pointers to index sets
3318 -  ov  - the additional overlap requested
3319 
3320    Level: developer
3321 
3322 .keywords: matrix, overlap, Schwarz
3323 
3324 .seealso: MatGetSubMatrices()
3325 @*/
3326 int MatIncreaseOverlap(Mat mat,int n,IS *is,int ov)
3327 {
3328   int ierr;
3329 
3330   PetscFunctionBegin;
3331   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3332   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3333   if (mat->factor)     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3334 
3335   if (!ov) PetscFunctionReturn(0);
3336   if (!mat->ops->increaseoverlap) SETERRQ(PETSC_ERR_SUP,0,"");
3337   PLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);
3338   ierr = (*mat->ops->increaseoverlap)(mat,n,is,ov);CHKERRQ(ierr);
3339   PLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);
3340   PetscFunctionReturn(0);
3341 }
3342 
3343 #undef __FUNC__
3344 #define __FUNC__ /*<a name=""></a>*/"MatPrintHelp"
3345 /*@
3346    MatPrintHelp - Prints all the options for the matrix.
3347 
3348    Collective on Mat
3349 
3350    Input Parameter:
3351 .  mat - the matrix
3352 
3353    Options Database Keys:
3354 +  -help - Prints matrix options
3355 -  -h - Prints matrix options
3356 
3357    Level: developer
3358 
3359 .keywords: mat, help
3360 
3361 .seealso: MatCreate(), MatCreateXXX()
3362 @*/
3363 int MatPrintHelp(Mat mat)
3364 {
3365   static PetscTruth called = PETSC_FALSE;
3366   int               ierr;
3367   MPI_Comm          comm;
3368 
3369   PetscFunctionBegin;
3370   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3371 
3372   comm = mat->comm;
3373   if (!called) {
3374     ierr = (*PetscHelpPrintf)(comm,"General matrix options:\n");CHKERRQ(ierr);
3375     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_info: view basic matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr);
3376     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_info_detailed: view detailed matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr);
3377     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_draw: draw nonzero matrix structure during MatAssemblyEnd()\n");CHKERRQ(ierr);
3378     ierr = (*PetscHelpPrintf)(comm,"      -draw_pause <sec>: set seconds of display pause\n");CHKERRQ(ierr);
3379     ierr = (*PetscHelpPrintf)(comm,"      -display <name>: set alternate display\n");CHKERRQ(ierr);
3380     called = PETSC_TRUE;
3381   }
3382   if (mat->ops->printhelp) {
3383     ierr = (*mat->ops->printhelp)(mat);CHKERRQ(ierr);
3384   }
3385   PetscFunctionReturn(0);
3386 }
3387 
3388 #undef __FUNC__
3389 #define __FUNC__ /*<a name=""></a>*/"MatGetBlockSize"
3390 /*@
3391    MatGetBlockSize - Returns the matrix block size; useful especially for the
3392    block row and block diagonal formats.
3393 
3394    Not Collective
3395 
3396    Input Parameter:
3397 .  mat - the matrix
3398 
3399    Output Parameter:
3400 .  bs - block size
3401 
3402    Notes:
3403    Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG.
3404    Block row formats are MATSEQBAIJ, MATMPIBAIJ
3405 
3406    Level: intermediate
3407 
3408 .keywords: matrix, get, block, size
3409 
3410 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag()
3411 @*/
3412 int MatGetBlockSize(Mat mat,int *bs)
3413 {
3414   int ierr;
3415 
3416   PetscFunctionBegin;
3417   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3418   PetscValidIntPointer(bs);
3419   if (!mat->ops->getblocksize) SETERRQ(PETSC_ERR_SUP,0,"");
3420   ierr = (*mat->ops->getblocksize)(mat,bs);CHKERRQ(ierr);
3421   PetscFunctionReturn(0);
3422 }
3423 
3424 #undef __FUNC__
3425 #define __FUNC__ /*<a name=""></a>*/"MatGetRowIJ"
3426 /*@C
3427     MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices.
3428 
3429    Collective on Mat
3430 
3431     Input Parameters:
3432 +   mat - the matrix
3433 .   shift -  0 or 1 indicating we want the indices starting at 0 or 1
3434 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3435                 symmetrized
3436 
3437     Output Parameters:
3438 +   n - number of rows in the (possibly compressed) matrix
3439 .   ia - the row pointers
3440 .   ja - the column indices
3441 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3442 
3443     Level: developer
3444 
3445 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3446 @*/
3447 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3448 {
3449   int ierr;
3450 
3451   PetscFunctionBegin;
3452   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3453   if (ia) PetscValidIntPointer(ia);
3454   if (ja) PetscValidIntPointer(ja);
3455   PetscValidIntPointer(done);
3456   if (!mat->ops->getrowij) *done = PETSC_FALSE;
3457   else {
3458     *done = PETSC_TRUE;
3459     ierr  = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3460   }
3461   PetscFunctionReturn(0);
3462 }
3463 
3464 #undef __FUNC__
3465 #define __FUNC__ /*<a name=""></a>*/"MatGetColumnIJ"
3466 /*@C
3467     MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices.
3468 
3469     Collective on Mat
3470 
3471     Input Parameters:
3472 +   mat - the matrix
3473 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3474 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3475                 symmetrized
3476 
3477     Output Parameters:
3478 +   n - number of columns in the (possibly compressed) matrix
3479 .   ia - the column pointers
3480 .   ja - the row indices
3481 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3482 
3483     Level: developer
3484 
3485 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3486 @*/
3487 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3488 {
3489   int ierr;
3490 
3491   PetscFunctionBegin;
3492   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3493   if (ia) PetscValidIntPointer(ia);
3494   if (ja) PetscValidIntPointer(ja);
3495   PetscValidIntPointer(done);
3496 
3497   if (!mat->ops->getcolumnij) *done = PETSC_FALSE;
3498   else {
3499     *done = PETSC_TRUE;
3500     ierr  = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3501   }
3502   PetscFunctionReturn(0);
3503 }
3504 
3505 #undef __FUNC__
3506 #define __FUNC__ /*<a name=""></a>*/"MatRestoreRowIJ"
3507 /*@C
3508     MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with
3509     MatGetRowIJ().
3510 
3511     Collective on Mat
3512 
3513     Input Parameters:
3514 +   mat - the matrix
3515 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3516 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3517                 symmetrized
3518 
3519     Output Parameters:
3520 +   n - size of (possibly compressed) matrix
3521 .   ia - the row pointers
3522 .   ja - the column indices
3523 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3524 
3525     Level: developer
3526 
3527 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3528 @*/
3529 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3530 {
3531   int ierr;
3532 
3533   PetscFunctionBegin;
3534   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3535   if (ia) PetscValidIntPointer(ia);
3536   if (ja) PetscValidIntPointer(ja);
3537   PetscValidIntPointer(done);
3538 
3539   if (!mat->ops->restorerowij) *done = PETSC_FALSE;
3540   else {
3541     *done = PETSC_TRUE;
3542     ierr  = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3543   }
3544   PetscFunctionReturn(0);
3545 }
3546 
3547 #undef __FUNC__
3548 #define __FUNC__ /*<a name=""></a>*/"MatRestoreColumnIJ"
3549 /*@C
3550     MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with
3551     MatGetColumnIJ().
3552 
3553     Collective on Mat
3554 
3555     Input Parameters:
3556 +   mat - the matrix
3557 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3558 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3559                 symmetrized
3560 
3561     Output Parameters:
3562 +   n - size of (possibly compressed) matrix
3563 .   ia - the column pointers
3564 .   ja - the row indices
3565 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3566 
3567     Level: developer
3568 
3569 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3570 @*/
3571 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3572 {
3573   int ierr;
3574 
3575   PetscFunctionBegin;
3576   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3577   if (ia) PetscValidIntPointer(ia);
3578   if (ja) PetscValidIntPointer(ja);
3579   PetscValidIntPointer(done);
3580 
3581   if (!mat->ops->restorecolumnij) *done = PETSC_FALSE;
3582   else {
3583     *done = PETSC_TRUE;
3584     ierr  = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3585   }
3586   PetscFunctionReturn(0);
3587 }
3588 
3589 #undef __FUNC__
3590 #define __FUNC__ /*<a name=""></a>*/"MatColoringPatch"
3591 /*@C
3592     MatColoringPatch -Used inside matrix coloring routines that
3593     use MatGetRowIJ() and/or MatGetColumnIJ().
3594 
3595     Collective on Mat
3596 
3597     Input Parameters:
3598 +   mat - the matrix
3599 .   n   - number of colors
3600 -   colorarray - array indicating color for each column
3601 
3602     Output Parameters:
3603 .   iscoloring - coloring generated using colorarray information
3604 
3605     Level: developer
3606 
3607 .seealso: MatGetRowIJ(), MatGetColumnIJ()
3608 
3609 .keywords: mat, coloring, patch
3610 @*/
3611 int MatColoringPatch(Mat mat,int n,int *colorarray,ISColoring *iscoloring)
3612 {
3613   int ierr;
3614 
3615   PetscFunctionBegin;
3616   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3617   PetscValidIntPointer(colorarray);
3618 
3619   if (!mat->ops->coloringpatch) {SETERRQ(PETSC_ERR_SUP,0,"");}
3620   else {
3621     ierr  = (*mat->ops->coloringpatch)(mat,n,colorarray,iscoloring);CHKERRQ(ierr);
3622   }
3623   PetscFunctionReturn(0);
3624 }
3625 
3626 
3627 #undef __FUNC__
3628 #define __FUNC__ /*<a name=""></a>*/"MatSetUnfactored"
3629 /*@
3630    MatSetUnfactored - Resets a factored matrix to be treated as unfactored.
3631 
3632    Collective on Mat
3633 
3634    Input Parameter:
3635 .  mat - the factored matrix to be reset
3636 
3637    Notes:
3638    This routine should be used only with factored matrices formed by in-place
3639    factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE
3640    format).  This option can save memory, for example, when solving nonlinear
3641    systems with a matrix-free Newton-Krylov method and a matrix-based, in-place
3642    ILU(0) preconditioner.
3643 
3644    Note that one can specify in-place ILU(0) factorization by calling
3645 .vb
3646      PCType(pc,PCILU);
3647      PCILUSeUseInPlace(pc);
3648 .ve
3649    or by using the options -pc_type ilu -pc_ilu_in_place
3650 
3651    In-place factorization ILU(0) can also be used as a local
3652    solver for the blocks within the block Jacobi or additive Schwarz
3653    methods (runtime option: -sub_pc_ilu_in_place).  See the discussion
3654    of these preconditioners in the users manual for details on setting
3655    local solver options.
3656 
3657    Most users should employ the simplified SLES interface for linear solvers
3658    instead of working directly with matrix algebra routines such as this.
3659    See, e.g., SLESCreate().
3660 
3661    Level: developer
3662 
3663 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace()
3664 
3665 .keywords: matrix-free, in-place ILU, in-place LU
3666 @*/
3667 int MatSetUnfactored(Mat mat)
3668 {
3669   int ierr;
3670 
3671   PetscFunctionBegin;
3672   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3673   mat->factor = 0;
3674   if (!mat->ops->setunfactored) PetscFunctionReturn(0);
3675   ierr = (*mat->ops->setunfactored)(mat);CHKERRQ(ierr);
3676   PetscFunctionReturn(0);
3677 }
3678 
3679 #undef __FUNC__
3680 #define __FUNC__ /*<a name=""></a>*/"MatGetType"
3681 /*@C
3682    MatGetType - Gets the matrix type and name (as a string) from the matrix.
3683 
3684    Not Collective
3685 
3686    Input Parameter:
3687 .  mat - the matrix
3688 
3689    Output Parameter:
3690 +  type - the matrix type (or use PETSC_NULL)
3691 -  name - name of matrix type (or use PETSC_NULL)
3692 
3693    Level: intermediate
3694 
3695 .keywords: matrix, get, type, name
3696 @*/
3697 int MatGetType(Mat mat,MatType *type,char **name)
3698 {
3699   int  itype = (int)mat->type;
3700   char *matname[10];
3701 
3702   PetscFunctionBegin;
3703   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3704 
3705   if (type) *type = (MatType) mat->type;
3706   if (name) {
3707     /* Note:  Be sure that this list corresponds to the enum in mat.h */
3708     matname[0] = "MATSEQDENSE";
3709     matname[1] = "MATSEQAIJ";
3710     matname[2] = "MATMPIAIJ";
3711     matname[3] = "MATSHELL";
3712     matname[4] = "MATMPIROWBS";
3713     matname[5] = "MATSEQBDIAG";
3714     matname[6] = "MATMPIBDIAG";
3715     matname[7] = "MATMPIDENSE";
3716     matname[8] = "MATSEQBAIJ";
3717     matname[9] = "MATMPIBAIJ";
3718 
3719     if (itype < 0 || itype > 9) *name = "Unknown matrix type";
3720     else                        *name = matname[itype];
3721   }
3722   PetscFunctionReturn(0);
3723 }
3724 
3725 /*MC
3726     MatGetArrayF90 - Accesses a matrix array from Fortran90.
3727 
3728     Synopsis:
3729     MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3730 
3731     Not collective
3732 
3733     Input Parameter:
3734 .   x - matrix
3735 
3736     Output Parameters:
3737 +   xx_v - the Fortran90 pointer to the array
3738 -   ierr - error code
3739 
3740     Example of Usage:
3741 .vb
3742       Scalar, pointer xx_v(:)
3743       ....
3744       call MatGetArrayF90(x,xx_v,ierr)
3745       a = xx_v(3)
3746       call MatRestoreArrayF90(x,xx_v,ierr)
3747 .ve
3748 
3749     Notes:
3750     Not yet supported for all F90 compilers
3751 
3752     Level: advanced
3753 
3754 .seealso:  MatRestoreArrayF90(), MatGetArray(), MatRestoreArray()
3755 
3756 .keywords:  matrix, array, f90
3757 M*/
3758 
3759 /*MC
3760     MatRestoreArrayF90 - Restores a matrix array that has been
3761     accessed with MatGetArrayF90().
3762 
3763     Synopsis:
3764     MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3765 
3766     Not collective
3767 
3768     Input Parameters:
3769 +   x - matrix
3770 -   xx_v - the Fortran90 pointer to the array
3771 
3772     Output Parameter:
3773 .   ierr - error code
3774 
3775     Example of Usage:
3776 .vb
3777        Scalar, pointer xx_v(:)
3778        ....
3779        call MatGetArrayF90(x,xx_v,ierr)
3780        a = xx_v(3)
3781        call MatRestoreArrayF90(x,xx_v,ierr)
3782 .ve
3783 
3784     Notes:
3785     Not yet supported for all F90 compilers
3786 
3787     Level: advanced
3788 
3789 .seealso:  MatGetArrayF90(), MatGetArray(), MatRestoreArray()
3790 
3791 .keywords:  matrix, array, f90
3792 M*/
3793 
3794 
3795 #undef __FUNC__
3796 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrix"
3797 /*@
3798     MatGetSubMatrix - Gets a single submatrix on the same number of processors
3799                       as the original matrix.
3800 
3801     Collective on Mat
3802 
3803     Input Parameters:
3804 +   mat - the original matrix
3805 .   isrow - rows this processor should obtain
3806 .   iscol - columns for all processors you wish to keep
3807 .   csize - number of columns "local" to this processor (does nothing for sequential
3808             matrices). This should match the result from VecGetLocalSize(x,...) if you
3809             plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE
3810 -   cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3811 
3812     Output Parameter:
3813 .   newmat - the new submatrix, of the same type as the old
3814 
3815     Level: advanced
3816 
3817 .keywords: matrix, get, submatrix, submatrices
3818 
3819 .seealso: MatGetSubMatrices()
3820 @*/
3821 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat)
3822 {
3823   int     ierr, size;
3824   Mat     *local;
3825 
3826   PetscFunctionBegin;
3827   ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr);
3828 
3829   /* if original matrix is on just one processor then use submatrix generated */
3830   if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) {
3831     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr);
3832     PetscFunctionReturn(0);
3833   } else if (!mat->ops->getsubmatrix && size == 1) {
3834     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr);
3835     *newmat = *local;
3836     ierr = PetscFree(local);CHKERRQ(ierr);
3837     PetscFunctionReturn(0);
3838   }
3839 
3840   if (!mat->ops->getsubmatrix) SETERRQ(PETSC_ERR_SUP,0,"Not currently implemented");
3841   ierr = (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr);
3842   PetscFunctionReturn(0);
3843 }
3844 
3845 #undef __FUNC__
3846 #define __FUNC__ /*<a name=""></a>*/"MatGetMaps"
3847 /*@C
3848    MatGetMaps - Returns the maps associated with the matrix.
3849 
3850    Not Collective
3851 
3852    Input Parameter:
3853 .  mat - the matrix
3854 
3855    Output Parameters:
3856 +  rmap - the row (right) map
3857 -  cmap - the column (left) map
3858 
3859    Level: developer
3860 
3861 .keywords: matrix, get, map
3862 @*/
3863 int MatGetMaps(Mat mat,Map *rmap,Map *cmap)
3864 {
3865   int ierr;
3866 
3867   PetscFunctionBegin;
3868   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3869   ierr = (*mat->ops->getmaps)(mat,rmap,cmap);CHKERRQ(ierr);
3870   PetscFunctionReturn(0);
3871 }
3872 
3873 /*
3874       Version that works for all PETSc matrices
3875 */
3876 #undef __FUNC__
3877 #define __FUNC__ /*<a name=""></a>*/"MatGetMaps_Petsc"
3878 int MatGetMaps_Petsc(Mat mat,Map *rmap,Map *cmap)
3879 {
3880   PetscFunctionBegin;
3881   if (rmap) *rmap = mat->rmap;
3882   if (cmap) *cmap = mat->cmap;
3883   PetscFunctionReturn(0);
3884 }
3885 
3886 #undef __FUNC__
3887 #define __FUNC__ /*<a name=""></a>*/"MatSetStashInitialSize"
3888 /*@
3889    MatSetStashInitialSize - sets the sizes of the matrix stash, that is
3890    used during the assembly process to store values that belong to
3891    other processors.
3892 
3893    Not Collective
3894 
3895    Input Parameters:
3896 +  mat   - the matrix
3897 .  size  - the initial size of the stash.
3898 -  bsize - the initial size of the block-stash(if used).
3899 
3900    Options Database Keys:
3901 +   -matstash_initial_size <size> or <size0,size1,...sizep-1>
3902 -   -matstash_block_initial_size <bsize>  or <bsize0,bsize1,...bsizep-1>
3903 
3904    Level: intermediate
3905 
3906    Notes:
3907      The block-stash is used for values set with VecSetValuesBlocked() while
3908      the stash is used for values set with VecSetValues()
3909 
3910      Run with the option -log_info and look for output of the form
3911      MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
3912      to determine the appropriate value, MM, to use for size and
3913      MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
3914      to determine the value, BMM to use for bsize
3915 
3916 .keywords: matrix, stash, assembly
3917 @*/
3918 int MatSetStashInitialSize(Mat mat,int size, int bsize)
3919 {
3920   int ierr;
3921 
3922   PetscFunctionBegin;
3923   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3924   ierr = MatStashSetInitialSize_Private(&mat->stash,size);CHKERRQ(ierr);
3925   ierr = MatStashSetInitialSize_Private(&mat->bstash,bsize);CHKERRQ(ierr);
3926   PetscFunctionReturn(0);
3927 }
3928 
3929 #undef __FUNC__
3930 #define __FUNC__ /*<a name=""></a>*/"MatInterpolateAdd"
3931 /*@
3932    MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of
3933      the matrix
3934 
3935    Collective on Mat
3936 
3937    Input Parameters:
3938 +  mat   - the matrix
3939 .  x,y - the vectors
3940 -  w - where the result is stored
3941 
3942    Level: intermediate
3943 
3944    Notes:
3945     w may be the same vector as y.
3946 
3947     This allows one to use either the restriction or interpolation (its transpose)
3948     matrix to do the interpolation
3949 
3950 .keywords: interpolate,
3951 
3952 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
3953 
3954 @*/
3955 int MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w)
3956 {
3957   int M,N,ierr;
3958 
3959   PetscFunctionBegin;
3960   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
3961   if (N > M) {
3962     ierr = MatMultTransposeAdd(A,x,y,w);CHKERRQ(ierr);
3963   } else {
3964     ierr = MatMultAdd(A,x,y,w);CHKERRQ(ierr);
3965   }
3966   PetscFunctionReturn(0);
3967 }
3968 
3969 #undef __FUNC__
3970 #define __FUNC__ /*<a name=""></a>*/"MatInterpolate"
3971 /*@
3972    MatInterpolate - y = A*x or A'*x depending on the shape of
3973      the matrix
3974 
3975    Collective on Mat
3976 
3977    Input Parameters:
3978 +  mat   - the matrix
3979 -  x,y - the vectors
3980 
3981    Level: intermediate
3982 
3983    Notes:
3984     This allows one to use either the restriction or interpolation (its transpose)
3985     matrix to do the interpolation
3986 
3987 .keywords: interpolate,
3988 
3989 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
3990 
3991 @*/
3992 int MatInterpolate(Mat A,Vec x,Vec y)
3993 {
3994   int M,N,ierr;
3995 
3996   PetscFunctionBegin;
3997   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
3998   if (N > M) {
3999     ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
4000   } else {
4001     ierr = MatMult(A,x,y);CHKERRQ(ierr);
4002   }
4003   PetscFunctionReturn(0);
4004 }
4005 
4006 #undef __FUNC__
4007 #define __FUNC__ /*<a name=""></a>*/"MatRestrict"
4008 /*@
4009    MatRestrict - y = A*x or A'*x
4010 
4011    Collective on Mat
4012 
4013    Input Parameters:
4014 +  mat   - the matrix
4015 -  x,y - the vectors
4016 
4017    Level: intermediate
4018 
4019    Notes:
4020     This allows one to use either the restriction or interpolation (its transpose)
4021     matrix to do the restriction
4022 
4023 .keywords: interpolate,
4024 
4025 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatInterpolate()
4026 
4027 @*/
4028 int MatRestrict(Mat A,Vec x,Vec y)
4029 {
4030   int M,N,ierr;
4031 
4032   PetscFunctionBegin;
4033   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
4034   if (N > M) {
4035     ierr = MatMult(A,x,y);CHKERRQ(ierr);
4036   } else {
4037     ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
4038   }
4039   PetscFunctionReturn(0);
4040 }
4041