xref: /petsc/src/mat/interface/matrix.c (revision 5ed6d96ad292ae16e6f6e2966d523a85f3a552b3)
1 /*$Id: matrix.c,v 1.363 2000/05/05 22:15:29 balay Exp bsmith $*/
2 
3 /*
4    This is where the abstract matrix operations are defined
5 */
6 
7 #include "src/mat/matimpl.h"        /*I "petscmat.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/petscmat.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/petscmat.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 -  fact - the matrix generated for the factor, from MatLUFactorSymbolic()
1335 
1336    Notes:
1337    See MatLUFactor() for in-place factorization.  See
1338    MatCholeskyFactorNumeric() for the symmetric, positive definite case.
1339 
1340    Most users should employ the simplified SLES interface for linear solvers
1341    instead of working directly with matrix algebra routines such as this.
1342    See, e.g., SLESCreate().
1343 
1344    Level: developer
1345 
1346 .keywords: matrix, factor, LU, numeric
1347 
1348 .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
1349 @*/
1350 int MatLUFactorNumeric(Mat mat,Mat *fact)
1351 {
1352   int        ierr;
1353   PetscTruth flg;
1354 
1355   PetscFunctionBegin;
1356   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1357   PetscValidPointer(fact);
1358   PetscValidHeaderSpecific(*fact,MAT_COOKIE);
1359   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1360   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1361     SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dimensions are different %d should = %d %d should = %d",
1362             mat->M,(*fact)->M,mat->N,(*fact)->N);
1363   }
1364   if (!(*fact)->ops->lufactornumeric) SETERRQ(PETSC_ERR_SUP,0,"");
1365 
1366   PLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0);
1367   ierr = (*(*fact)->ops->lufactornumeric)(mat,fact);CHKERRQ(ierr);
1368   PLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0);
1369   ierr = OptionsHasName(PETSC_NULL,"-mat_view_draw",&flg);CHKERRQ(ierr);
1370   if (flg) {
1371     ierr = OptionsHasName(PETSC_NULL,"-mat_view_contour",&flg);CHKERRQ(ierr);
1372     if (flg) {
1373       ViewerPushFormat(VIEWER_DRAW_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr);
1374     }
1375     ierr = MatView(*fact,VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
1376     ierr = ViewerFlush(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
1377     if (flg) {
1378       ViewerPopFormat(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
1379     }
1380   }
1381   PetscFunctionReturn(0);
1382 }
1383 
1384 #undef __FUNC__
1385 #define __FUNC__ /*<a name=""></a>*/"MatCholeskyFactor"
1386 /*@
1387    MatCholeskyFactor - Performs in-place Cholesky factorization of a
1388    symmetric matrix.
1389 
1390    Collective on Mat
1391 
1392    Input Parameters:
1393 +  mat - the matrix
1394 .  perm - row and column permutations
1395 -  f - expected fill as ratio of original fill
1396 
1397    Notes:
1398    See MatLUFactor() for the nonsymmetric case.  See also
1399    MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().
1400 
1401    Most users should employ the simplified SLES interface for linear solvers
1402    instead of working directly with matrix algebra routines such as this.
1403    See, e.g., SLESCreate().
1404 
1405    Level: developer
1406 
1407 .keywords: matrix, factor, in-place, Cholesky
1408 
1409 .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
1410           MatGetOrdering()
1411 
1412 @*/
1413 int MatCholeskyFactor(Mat mat,IS perm,PetscReal f)
1414 {
1415   int ierr;
1416 
1417   PetscFunctionBegin;
1418   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1419   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Matrix must be square");
1420   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1421   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1422   if (!mat->ops->choleskyfactor) SETERRQ(PETSC_ERR_SUP,0,"");
1423 
1424   PLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);
1425   ierr = (*mat->ops->choleskyfactor)(mat,perm,f);CHKERRQ(ierr);
1426   PLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);
1427   PetscFunctionReturn(0);
1428 }
1429 
1430 #undef __FUNC__
1431 #define __FUNC__ /*<a name=""></a>*/"MatCholeskyFactorSymbolic"
1432 /*@
1433    MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
1434    of a symmetric matrix.
1435 
1436    Collective on Mat
1437 
1438    Input Parameters:
1439 +  mat - the matrix
1440 .  perm - row and column permutations
1441 -  f - expected fill as ratio of original
1442 
1443    Output Parameter:
1444 .  fact - the factored matrix
1445 
1446    Notes:
1447    See MatLUFactorSymbolic() for the nonsymmetric case.  See also
1448    MatCholeskyFactor() and MatCholeskyFactorNumeric().
1449 
1450    Most users should employ the simplified SLES interface for linear solvers
1451    instead of working directly with matrix algebra routines such as this.
1452    See, e.g., SLESCreate().
1453 
1454    Level: developer
1455 
1456 .keywords: matrix, factor, factorization, symbolic, Cholesky
1457 
1458 .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
1459           MatGetOrdering()
1460 
1461 @*/
1462 int MatCholeskyFactorSymbolic(Mat mat,IS perm,PetscReal f,Mat *fact)
1463 {
1464   int ierr;
1465 
1466   PetscFunctionBegin;
1467   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1468   PetscValidPointer(fact);
1469   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Matrix must be square");
1470   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1471   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1472   if (!mat->ops->choleskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"");
1473 
1474   PLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
1475   ierr = (*mat->ops->choleskyfactorsymbolic)(mat,perm,f,fact);CHKERRQ(ierr);
1476   PLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
1477   PetscFunctionReturn(0);
1478 }
1479 
1480 #undef __FUNC__
1481 #define __FUNC__ /*<a name=""></a>*/"MatCholeskyFactorNumeric"
1482 /*@
1483    MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
1484    of a symmetric matrix. Call this routine after first calling
1485    MatCholeskyFactorSymbolic().
1486 
1487    Collective on Mat
1488 
1489    Input Parameter:
1490 .  mat - the initial matrix
1491 
1492    Output Parameter:
1493 .  fact - the factored matrix
1494 
1495    Notes:
1496    Most users should employ the simplified SLES interface for linear solvers
1497    instead of working directly with matrix algebra routines such as this.
1498    See, e.g., SLESCreate().
1499 
1500    Level: developer
1501 
1502 .keywords: matrix, factor, numeric, Cholesky
1503 
1504 .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
1505 @*/
1506 int MatCholeskyFactorNumeric(Mat mat,Mat *fact)
1507 {
1508   int ierr;
1509 
1510   PetscFunctionBegin;
1511   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1512   PetscValidPointer(fact);
1513   if (!mat->ops->choleskyfactornumeric) SETERRQ(PETSC_ERR_SUP,0,"");
1514   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1515   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1516     SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dim %d should = %d %d should = %d",
1517             mat->M,(*fact)->M,mat->N,(*fact)->N);
1518   }
1519 
1520   PLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
1521   ierr = (*mat->ops->choleskyfactornumeric)(mat,fact);CHKERRQ(ierr);
1522   PLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
1523   PetscFunctionReturn(0);
1524 }
1525 
1526 /* ----------------------------------------------------------------*/
1527 #undef __FUNC__
1528 #define __FUNC__ /*<a name=""></a>*/"MatSolve"
1529 /*@
1530    MatSolve - Solves A x = b, given a factored matrix.
1531 
1532    Collective on Mat and Vec
1533 
1534    Input Parameters:
1535 +  mat - the factored matrix
1536 -  b - the right-hand-side vector
1537 
1538    Output Parameter:
1539 .  x - the result vector
1540 
1541    Notes:
1542    The vectors b and x cannot be the same.  I.e., one cannot
1543    call MatSolve(A,x,x).
1544 
1545    Notes:
1546    Most users should employ the simplified SLES interface for linear solvers
1547    instead of working directly with matrix algebra routines such as this.
1548    See, e.g., SLESCreate().
1549 
1550    Level: developer
1551 
1552 .keywords: matrix, linear system, solve, LU, Cholesky, triangular solve
1553 
1554 .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd()
1555 @*/
1556 int MatSolve(Mat mat,Vec b,Vec x)
1557 {
1558   int ierr;
1559 
1560   PetscFunctionBegin;
1561   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1562   PetscValidHeaderSpecific(b,VEC_COOKIE);
1563   PetscValidHeaderSpecific(x,VEC_COOKIE);
1564   PetscCheckSameComm(mat,b);
1565   PetscCheckSameComm(mat,x);
1566   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1567   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1568   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1569   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1570   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1571   if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0);
1572 
1573   if (!mat->ops->solve) SETERRQ(PETSC_ERR_SUP,0,"");
1574   PLogEventBegin(MAT_Solve,mat,b,x,0);
1575   ierr = (*mat->ops->solve)(mat,b,x);CHKERRQ(ierr);
1576   PLogEventEnd(MAT_Solve,mat,b,x,0);
1577   PetscFunctionReturn(0);
1578 }
1579 
1580 #undef __FUNC__
1581 #define __FUNC__ /*<a name=""></a>*/"MatForwardSolve"
1582 /* @
1583    MatForwardSolve - Solves L x = b, given a factored matrix, A = LU.
1584 
1585    Collective on Mat and Vec
1586 
1587    Input Parameters:
1588 +  mat - the factored matrix
1589 -  b - the right-hand-side vector
1590 
1591    Output Parameter:
1592 .  x - the result vector
1593 
1594    Notes:
1595    MatSolve() should be used for most applications, as it performs
1596    a forward solve followed by a backward solve.
1597 
1598    The vectors b and x cannot be the same.  I.e., one cannot
1599    call MatForwardSolve(A,x,x).
1600 
1601    Most users should employ the simplified SLES interface for linear solvers
1602    instead of working directly with matrix algebra routines such as this.
1603    See, e.g., SLESCreate().
1604 
1605    Level: developer
1606 
1607 .keywords: matrix, forward, LU, Cholesky, triangular solve
1608 
1609 .seealso: MatSolve(), MatBackwardSolve()
1610 @ */
1611 int MatForwardSolve(Mat mat,Vec b,Vec x)
1612 {
1613   int ierr;
1614 
1615   PetscFunctionBegin;
1616   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1617   PetscValidHeaderSpecific(b,VEC_COOKIE);
1618   PetscValidHeaderSpecific(x,VEC_COOKIE);
1619   PetscCheckSameComm(mat,b);
1620   PetscCheckSameComm(mat,x);
1621   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1622   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1623   if (!mat->ops->forwardsolve) SETERRQ(PETSC_ERR_SUP,0,"");
1624   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1625   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1626   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1627 
1628   PLogEventBegin(MAT_ForwardSolve,mat,b,x,0);
1629   ierr = (*mat->ops->forwardsolve)(mat,b,x);CHKERRQ(ierr);
1630   PLogEventEnd(MAT_ForwardSolve,mat,b,x,0);
1631   PetscFunctionReturn(0);
1632 }
1633 
1634 #undef __FUNC__
1635 #define __FUNC__ /*<a name=""></a>*/"MatBackwardSolve"
1636 /* @
1637    MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.
1638 
1639    Collective on Mat and Vec
1640 
1641    Input Parameters:
1642 +  mat - the factored matrix
1643 -  b - the right-hand-side vector
1644 
1645    Output Parameter:
1646 .  x - the result vector
1647 
1648    Notes:
1649    MatSolve() should be used for most applications, as it performs
1650    a forward solve followed by a backward solve.
1651 
1652    The vectors b and x cannot be the same.  I.e., one cannot
1653    call MatBackwardSolve(A,x,x).
1654 
1655    Most users should employ the simplified SLES interface for linear solvers
1656    instead of working directly with matrix algebra routines such as this.
1657    See, e.g., SLESCreate().
1658 
1659    Level: developer
1660 
1661 .keywords: matrix, backward, LU, Cholesky, triangular solve
1662 
1663 .seealso: MatSolve(), MatForwardSolve()
1664 @ */
1665 int MatBackwardSolve(Mat mat,Vec b,Vec x)
1666 {
1667   int ierr;
1668 
1669   PetscFunctionBegin;
1670   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1671   PetscValidHeaderSpecific(b,VEC_COOKIE);
1672   PetscValidHeaderSpecific(x,VEC_COOKIE);
1673   PetscCheckSameComm(mat,b);
1674   PetscCheckSameComm(mat,x);
1675   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1676   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1677   if (!mat->ops->backwardsolve) SETERRQ(PETSC_ERR_SUP,0,"");
1678   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1679   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1680   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1681 
1682   PLogEventBegin(MAT_BackwardSolve,mat,b,x,0);
1683   ierr = (*mat->ops->backwardsolve)(mat,b,x);CHKERRQ(ierr);
1684   PLogEventEnd(MAT_BackwardSolve,mat,b,x,0);
1685   PetscFunctionReturn(0);
1686 }
1687 
1688 #undef __FUNC__
1689 #define __FUNC__ /*<a name=""></a>*/"MatSolveAdd"
1690 /*@
1691    MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.
1692 
1693    Collective on Mat and Vec
1694 
1695    Input Parameters:
1696 +  mat - the factored matrix
1697 .  b - the right-hand-side vector
1698 -  y - the vector to be added to
1699 
1700    Output Parameter:
1701 .  x - the result vector
1702 
1703    Notes:
1704    The vectors b and x cannot be the same.  I.e., one cannot
1705    call MatSolveAdd(A,x,y,x).
1706 
1707    Most users should employ the simplified SLES interface for linear solvers
1708    instead of working directly with matrix algebra routines such as this.
1709    See, e.g., SLESCreate().
1710 
1711    Level: developer
1712 
1713 .keywords: matrix, linear system, solve, LU, Cholesky, add
1714 
1715 .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd()
1716 @*/
1717 int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
1718 {
1719   Scalar one = 1.0;
1720   Vec    tmp;
1721   int    ierr;
1722 
1723   PetscFunctionBegin;
1724   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1725   PetscValidHeaderSpecific(y,VEC_COOKIE);
1726   PetscValidHeaderSpecific(b,VEC_COOKIE);
1727   PetscValidHeaderSpecific(x,VEC_COOKIE);
1728   PetscCheckSameComm(mat,b);
1729   PetscCheckSameComm(mat,y);
1730   PetscCheckSameComm(mat,x);
1731   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1732   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1733   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1734   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1735   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
1736   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1737   if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim %d %d",x->n,y->n);
1738 
1739   PLogEventBegin(MAT_SolveAdd,mat,b,x,y);
1740   if (mat->ops->solveadd)  {
1741     ierr = (*mat->ops->solveadd)(mat,b,y,x);CHKERRQ(ierr);
1742   } else {
1743     /* do the solve then the add manually */
1744     if (x != y) {
1745       ierr = MatSolve(mat,b,x);CHKERRQ(ierr);
1746       ierr = VecAXPY(&one,y,x);CHKERRQ(ierr);
1747     } else {
1748       ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr);
1749       PLogObjectParent(mat,tmp);
1750       ierr = VecCopy(x,tmp);CHKERRQ(ierr);
1751       ierr = MatSolve(mat,b,x);CHKERRQ(ierr);
1752       ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr);
1753       ierr = VecDestroy(tmp);CHKERRQ(ierr);
1754     }
1755   }
1756   PLogEventEnd(MAT_SolveAdd,mat,b,x,y);
1757   PetscFunctionReturn(0);
1758 }
1759 
1760 #undef __FUNC__
1761 #define __FUNC__ /*<a name=""></a>*/"MatSolveTranspose"
1762 /*@
1763    MatSolveTranspose - Solves A' x = b, given a factored matrix.
1764 
1765    Collective on Mat and Vec
1766 
1767    Input Parameters:
1768 +  mat - the factored matrix
1769 -  b - the right-hand-side vector
1770 
1771    Output Parameter:
1772 .  x - the result vector
1773 
1774    Notes:
1775    The vectors b and x cannot be the same.  I.e., one cannot
1776    call MatSolveTranspose(A,x,x).
1777 
1778    Most users should employ the simplified SLES interface for linear solvers
1779    instead of working directly with matrix algebra routines such as this.
1780    See, e.g., SLESCreate().
1781 
1782    Level: developer
1783 
1784 .keywords: matrix, linear system, solve, LU, Cholesky, transpose
1785 
1786 .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd()
1787 @*/
1788 int MatSolveTranspose(Mat mat,Vec b,Vec x)
1789 {
1790   int ierr;
1791 
1792   PetscFunctionBegin;
1793   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1794   PetscValidHeaderSpecific(b,VEC_COOKIE);
1795   PetscValidHeaderSpecific(x,VEC_COOKIE);
1796   PetscCheckSameComm(mat,b);
1797   PetscCheckSameComm(mat,x);
1798   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1799   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1800   if (!mat->ops->solvetranspose) SETERRQ(PETSC_ERR_SUP,0,"");
1801   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
1802   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);
1803 
1804   PLogEventBegin(MAT_SolveTranspose,mat,b,x,0);
1805   ierr = (*mat->ops->solvetranspose)(mat,b,x);CHKERRQ(ierr);
1806   PLogEventEnd(MAT_SolveTranspose,mat,b,x,0);
1807   PetscFunctionReturn(0);
1808 }
1809 
1810 #undef __FUNC__
1811 #define __FUNC__ /*<a name=""></a>*/"MatSolveTransposeAdd"
1812 /*@
1813    MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a
1814                       factored matrix.
1815 
1816    Collective on Mat and Vec
1817 
1818    Input Parameters:
1819 +  mat - the factored matrix
1820 .  b - the right-hand-side vector
1821 -  y - the vector to be added to
1822 
1823    Output Parameter:
1824 .  x - the result vector
1825 
1826    Notes:
1827    The vectors b and x cannot be the same.  I.e., one cannot
1828    call MatSolveTransposeAdd(A,x,y,x).
1829 
1830    Most users should employ the simplified SLES interface for linear solvers
1831    instead of working directly with matrix algebra routines such as this.
1832    See, e.g., SLESCreate().
1833 
1834    Level: developer
1835 
1836 .keywords: matrix, linear system, solve, LU, Cholesky, transpose, add
1837 
1838 .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose()
1839 @*/
1840 int MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x)
1841 {
1842   Scalar one = 1.0;
1843   int    ierr;
1844   Vec    tmp;
1845 
1846   PetscFunctionBegin;
1847   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1848   PetscValidHeaderSpecific(y,VEC_COOKIE);
1849   PetscValidHeaderSpecific(b,VEC_COOKIE);
1850   PetscValidHeaderSpecific(x,VEC_COOKIE);
1851   PetscCheckSameComm(mat,b);
1852   PetscCheckSameComm(mat,y);
1853   PetscCheckSameComm(mat,x);
1854   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1855   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1856   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
1857   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);
1858   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
1859   if (x->n != y->n)   SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim %d %d",x->n,y->n);
1860 
1861   PLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);
1862   if (mat->ops->solvetransposeadd) {
1863     ierr = (*mat->ops->solvetransposeadd)(mat,b,y,x);CHKERRQ(ierr);
1864   } else {
1865     /* do the solve then the add manually */
1866     if (x != y) {
1867       ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr);
1868       ierr = VecAXPY(&one,y,x);CHKERRQ(ierr);
1869     } else {
1870       ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr);
1871       PLogObjectParent(mat,tmp);
1872       ierr = VecCopy(x,tmp);CHKERRQ(ierr);
1873       ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr);
1874       ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr);
1875       ierr = VecDestroy(tmp);CHKERRQ(ierr);
1876     }
1877   }
1878   PLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);
1879   PetscFunctionReturn(0);
1880 }
1881 /* ----------------------------------------------------------------*/
1882 
1883 #undef __FUNC__
1884 #define __FUNC__ /*<a name=""></a>*/"MatRelax"
1885 /*@
1886    MatRelax - Computes one relaxation sweep.
1887 
1888    Collective on Mat and Vec
1889 
1890    Input Parameters:
1891 +  mat - the matrix
1892 .  b - the right hand side
1893 .  omega - the relaxation factor
1894 .  flag - flag indicating the type of SOR (see below)
1895 .  shift -  diagonal shift
1896 -  its - the number of iterations
1897 
1898    Output Parameters:
1899 .  x - the solution (can contain an initial guess)
1900 
1901    SOR Flags:
1902 .     SOR_FORWARD_SWEEP - forward SOR
1903 .     SOR_BACKWARD_SWEEP - backward SOR
1904 .     SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
1905 .     SOR_LOCAL_FORWARD_SWEEP - local forward SOR
1906 .     SOR_LOCAL_BACKWARD_SWEEP - local forward SOR
1907 .     SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
1908 .     SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies
1909          upper/lower triangular part of matrix to
1910          vector (with omega)
1911 .     SOR_ZERO_INITIAL_GUESS - zero initial guess
1912 
1913    Notes:
1914    SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
1915    SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings
1916    on each processor.
1917 
1918    Application programmers will not generally use MatRelax() directly,
1919    but instead will employ the SLES/PC interface.
1920 
1921    Notes for Advanced Users:
1922    The flags are implemented as bitwise inclusive or operations.
1923    For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
1924    to specify a zero initial guess for SSOR.
1925 
1926    Most users should employ the simplified SLES interface for linear solvers
1927    instead of working directly with matrix algebra routines such as this.
1928    See, e.g., SLESCreate().
1929 
1930    Level: developer
1931 
1932 .keywords: matrix, relax, relaxation, sweep
1933 @*/
1934 int MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,int its,Vec x)
1935 {
1936   int ierr;
1937 
1938   PetscFunctionBegin;
1939   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1940   PetscValidHeaderSpecific(b,VEC_COOKIE);
1941   PetscValidHeaderSpecific(x,VEC_COOKIE);
1942   PetscCheckSameComm(mat,b);
1943   PetscCheckSameComm(mat,x);
1944   if (!mat->ops->relax) SETERRQ(PETSC_ERR_SUP,0,"");
1945   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1946   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1947   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1948   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1949   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1950 
1951   PLogEventBegin(MAT_Relax,mat,b,x,0);
1952   ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,x);CHKERRQ(ierr);
1953   PLogEventEnd(MAT_Relax,mat,b,x,0);
1954   PetscFunctionReturn(0);
1955 }
1956 
1957 #undef __FUNC__
1958 #define __FUNC__ /*<a name=""></a>*/"MatCopy_Basic"
1959 /*
1960       Default matrix copy routine.
1961 */
1962 int MatCopy_Basic(Mat A,Mat B,MatStructure str)
1963 {
1964   int    ierr,i,rstart,rend,nz,*cwork;
1965   Scalar *vwork;
1966 
1967   PetscFunctionBegin;
1968   ierr = MatZeroEntries(B);CHKERRQ(ierr);
1969   ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
1970   for (i=rstart; i<rend; i++) {
1971     ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr);
1972     ierr = MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);CHKERRQ(ierr);
1973     ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr);
1974   }
1975   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1976   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1977   PetscFunctionReturn(0);
1978 }
1979 
1980 #undef __FUNC__
1981 #define __FUNC__ /*<a name=""></a>*/"MatCopy"
1982 /*@C
1983    MatCopy - Copys a matrix to another matrix.
1984 
1985    Collective on Mat
1986 
1987    Input Parameters:
1988 +  A - the matrix
1989 -  str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN
1990 
1991    Output Parameter:
1992 .  B - where the copy is put
1993 
1994    Notes:
1995    If you use SAME_NONZERO_PATTERN then the zero matrices had better have the
1996    same nonzero pattern or the routine will crash.
1997 
1998    MatCopy() copies the matrix entries of a matrix to another existing
1999    matrix (after first zeroing the second matrix).  A related routine is
2000    MatConvert(), which first creates a new matrix and then copies the data.
2001 
2002    Level: intermediate
2003 
2004 .keywords: matrix, copy, convert
2005 
2006 .seealso: MatConvert()
2007 @*/
2008 int MatCopy(Mat A,Mat B,MatStructure str)
2009 {
2010   int ierr;
2011 
2012   PetscFunctionBegin;
2013   PetscValidHeaderSpecific(A,MAT_COOKIE);
2014   PetscValidHeaderSpecific(B,MAT_COOKIE);
2015   PetscCheckSameComm(A,B);
2016   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2017   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2018   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,
2019                                              A->N,B->N);
2020 
2021   PLogEventBegin(MAT_Copy,A,B,0,0);
2022   if (A->ops->copy) {
2023     ierr = (*A->ops->copy)(A,B,str);CHKERRQ(ierr);
2024   } else { /* generic conversion */
2025     ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr);
2026   }
2027   PLogEventEnd(MAT_Copy,A,B,0,0);
2028   PetscFunctionReturn(0);
2029 }
2030 
2031 static int MatConvertersSet = 0;
2032 static int (*MatConverters[MAX_MATRIX_TYPES][MAX_MATRIX_TYPES])(Mat,MatType,Mat*) =
2033            {{0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
2034             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
2035             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
2036             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
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 
2040 #undef __FUNC__
2041 #define __FUNC__ /*<a name=""></a>*/"MatConvertRegister"
2042 /*@C
2043     MatConvertRegister - Allows one to register a routine that converts between
2044     two matrix types.
2045 
2046     Not Collective
2047 
2048     Input Parameters:
2049 +   intype - the type of matrix (defined in include/petscmat.h), for example, MATSEQAIJ.
2050 -   outtype - new matrix type, or MATSAME
2051 
2052     Level: advanced
2053 
2054 .seealso: MatConvertRegisterAll()
2055 @*/
2056 int MatConvertRegister(MatType intype,MatType outtype,int (*converter)(Mat,MatType,Mat*))
2057 {
2058   PetscFunctionBegin;
2059   MatConverters[intype][outtype] = converter;
2060   MatConvertersSet               = 1;
2061   PetscFunctionReturn(0);
2062 }
2063 
2064 #undef __FUNC__
2065 #define __FUNC__ /*<a name=""></a>*/"MatConvert"
2066 /*@C
2067    MatConvert - Converts a matrix to another matrix, either of the same
2068    or different type.
2069 
2070    Collective on Mat
2071 
2072    Input Parameters:
2073 +  mat - the matrix
2074 -  newtype - new matrix type.  Use MATSAME to create a new matrix of the
2075    same type as the original matrix.
2076 
2077    Output Parameter:
2078 .  M - pointer to place new matrix
2079 
2080    Notes:
2081    MatConvert() first creates a new matrix and then copies the data from
2082    the first matrix.  A related routine is MatCopy(), which copies the matrix
2083    entries of one matrix to another already existing matrix context.
2084 
2085    Level: intermediate
2086 
2087 .keywords: matrix, copy, convert
2088 
2089 .seealso: MatCopy(), MatDuplicate()
2090 @*/
2091 int MatConvert(Mat mat,MatType newtype,Mat *M)
2092 {
2093   int ierr;
2094 
2095   PetscFunctionBegin;
2096   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2097   PetscValidPointer(M);
2098   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2099   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2100 
2101   if (newtype > MAX_MATRIX_TYPES || newtype < -1) {
2102     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,1,"Not a valid matrix type");
2103   }
2104   *M  = 0;
2105 
2106   if (!MatConvertersSet) {
2107     ierr = MatLoadRegisterAll();CHKERRQ(ierr);
2108   }
2109 
2110   PLogEventBegin(MAT_Convert,mat,0,0,0);
2111   if ((newtype == mat->type || newtype == MATSAME) && mat->ops->duplicate) {
2112     ierr = (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);CHKERRQ(ierr);
2113   } else {
2114     if (!MatConvertersSet) {
2115       ierr = MatConvertRegisterAll();CHKERRQ(ierr);
2116     }
2117     if (!MatConverters[mat->type][newtype]) {
2118       SETERRQ(PETSC_ERR_ARG_WRONG,1,"Invalid matrix type, or matrix converter not registered");
2119     }
2120     ierr = (*MatConverters[mat->type][newtype])(mat,newtype,M);CHKERRQ(ierr);
2121   }
2122   PLogEventEnd(MAT_Convert,mat,0,0,0);
2123   PetscFunctionReturn(0);
2124 }
2125 
2126 #undef __FUNC__
2127 #define __FUNC__ /*<a name=""></a>*/"MatDuplicate"
2128 /*@C
2129    MatDuplicate - Duplicates a matrix including the non-zero structure.
2130 
2131    Collective on Mat
2132 
2133    Input Parameters:
2134 +  mat - the matrix
2135 -  op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero
2136         values as well or not
2137 
2138    Output Parameter:
2139 .  M - pointer to place new matrix
2140 
2141    Level: intermediate
2142 
2143 .keywords: matrix, copy, convert, duplicate
2144 
2145 .seealso: MatCopy(), MatConvert()
2146 @*/
2147 int MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
2148 {
2149   int ierr;
2150 
2151   PetscFunctionBegin;
2152   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2153   PetscValidPointer(M);
2154   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2155   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2156 
2157   *M  = 0;
2158   PLogEventBegin(MAT_Convert,mat,0,0,0);
2159   if (!mat->ops->duplicate) {
2160     SETERRQ(PETSC_ERR_SUP,1,"Not written for this matrix type");
2161   }
2162   ierr = (*mat->ops->duplicate)(mat,op,M);CHKERRQ(ierr);
2163   PLogEventEnd(MAT_Convert,mat,0,0,0);
2164   PetscFunctionReturn(0);
2165 }
2166 
2167 #undef __FUNC__
2168 #define __FUNC__ /*<a name=""></a>*/"MatGetDiagonal"
2169 /*@
2170    MatGetDiagonal - Gets the diagonal of a matrix.
2171 
2172    Collective on Mat and Vec
2173 
2174    Input Parameters:
2175 +  mat - the matrix
2176 -  v - the vector for storing the diagonal
2177 
2178    Output Parameter:
2179 .  v - the diagonal of the matrix
2180 
2181    Notes:
2182    For the SeqAIJ matrix format, this routine may also be called
2183    on a LU factored matrix; in that case it routines the reciprocal of
2184    the diagonal entries in U. It returns the entries permuted by the
2185    row and column permutation used during the symbolic factorization.
2186 
2187    Level: intermediate
2188 
2189 .keywords: matrix, get, diagonal
2190 
2191 .seealso: MatGetRow(), MatGetSubmatrices(), MatGetSubmatrix()
2192 @*/
2193 int MatGetDiagonal(Mat mat,Vec v)
2194 {
2195   int ierr;
2196 
2197   PetscFunctionBegin;
2198   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2199   PetscValidHeaderSpecific(v,VEC_COOKIE);
2200   /* PetscCheckSameComm(mat,v); Could be MPI vector but Seq matrix cause of two submatrix storage */
2201   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2202   if (!mat->ops->getdiagonal) SETERRQ(PETSC_ERR_SUP,0,"");
2203 
2204   ierr = (*mat->ops->getdiagonal)(mat,v);CHKERRQ(ierr);
2205   PetscFunctionReturn(0);
2206 }
2207 
2208 #undef __FUNC__
2209 #define __FUNC__ /*<a name=""></a>*/"MatTranspose"
2210 /*@C
2211    MatTranspose - Computes an in-place or out-of-place transpose of a matrix.
2212 
2213    Collective on Mat
2214 
2215    Input Parameter:
2216 .  mat - the matrix to transpose
2217 
2218    Output Parameters:
2219 .  B - the transpose (or pass in PETSC_NULL for an in-place transpose)
2220 
2221    Level: intermediate
2222 
2223 .keywords: matrix, transpose
2224 
2225 .seealso: MatMultTranspose(), MatMultTransposeAdd()
2226 @*/
2227 int MatTranspose(Mat mat,Mat *B)
2228 {
2229   int ierr;
2230 
2231   PetscFunctionBegin;
2232   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2233   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2234   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2235   if (!mat->ops->transpose) SETERRQ(PETSC_ERR_SUP,0,"");
2236   ierr = (*mat->ops->transpose)(mat,B);CHKERRQ(ierr);
2237   PetscFunctionReturn(0);
2238 }
2239 
2240 #undef __FUNC__
2241 #define __FUNC__ /*<a name=""></a>*/"MatPermute"
2242 /*@C
2243    MatPermute - Creates a new matrix with rows and columns permuted from the
2244    original.
2245 
2246    Collective on Mat
2247 
2248    Input Parameters:
2249 +  mat - the matrix to permute
2250 .  row - row permutation, each processor supplies only the permutation for its rows
2251 -  col - column permutation, each processor needs the entire column permutation, that is
2252          this is the same size as the total number of columns in the matrix
2253 
2254    Output Parameters:
2255 .  B - the permuted matrix
2256 
2257    Level: advanced
2258 
2259 .keywords: matrix, transpose
2260 
2261 .seealso: MatGetOrdering()
2262 @*/
2263 int MatPermute(Mat mat,IS row,IS col,Mat *B)
2264 {
2265   int ierr;
2266 
2267   PetscFunctionBegin;
2268   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2269   PetscValidHeaderSpecific(row,IS_COOKIE);
2270   PetscValidHeaderSpecific(col,IS_COOKIE);
2271   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2272   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2273   if (!mat->ops->permute) SETERRQ(PETSC_ERR_SUP,0,"");
2274   ierr = (*mat->ops->permute)(mat,row,col,B);CHKERRQ(ierr);
2275   PetscFunctionReturn(0);
2276 }
2277 
2278 #undef __FUNC__
2279 #define __FUNC__ /*<a name=""></a>*/"MatEqual"
2280 /*@
2281    MatEqual - Compares two matrices.
2282 
2283    Collective on Mat
2284 
2285    Input Parameters:
2286 +  A - the first matrix
2287 -  B - the second matrix
2288 
2289    Output Parameter:
2290 .  flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.
2291 
2292    Level: intermediate
2293 
2294 .keywords: matrix, equal, equivalent
2295 @*/
2296 int MatEqual(Mat A,Mat B,PetscTruth *flg)
2297 {
2298   int ierr;
2299 
2300   PetscFunctionBegin;
2301   PetscValidHeaderSpecific(A,MAT_COOKIE);
2302   PetscValidHeaderSpecific(B,MAT_COOKIE);
2303   PetscValidIntPointer(flg);
2304   PetscCheckSameComm(A,B);
2305   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2306   if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2307   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",
2308                                              A->M,B->M,A->N,B->N);
2309   if (!A->ops->equal) SETERRQ(PETSC_ERR_SUP,0,"");
2310   ierr = (*A->ops->equal)(A,B,flg);CHKERRQ(ierr);
2311   PetscFunctionReturn(0);
2312 }
2313 
2314 #undef __FUNC__
2315 #define __FUNC__ /*<a name=""></a>*/"MatDiagonalScale"
2316 /*@
2317    MatDiagonalScale - Scales a matrix on the left and right by diagonal
2318    matrices that are stored as vectors.  Either of the two scaling
2319    matrices can be PETSC_NULL.
2320 
2321    Collective on Mat
2322 
2323    Input Parameters:
2324 +  mat - the matrix to be scaled
2325 .  l - the left scaling vector (or PETSC_NULL)
2326 -  r - the right scaling vector (or PETSC_NULL)
2327 
2328    Notes:
2329    MatDiagonalScale() computes A = LAR, where
2330    L = a diagonal matrix, R = a diagonal matrix
2331 
2332    Level: intermediate
2333 
2334 .keywords: matrix, diagonal, scale
2335 
2336 .seealso: MatScale()
2337 @*/
2338 int MatDiagonalScale(Mat mat,Vec l,Vec r)
2339 {
2340   int ierr;
2341 
2342   PetscFunctionBegin;
2343   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2344   if (!mat->ops->diagonalscale) SETERRQ(PETSC_ERR_SUP,0,"");
2345   if (l) {PetscValidHeaderSpecific(l,VEC_COOKIE);PetscCheckSameComm(mat,l);}
2346   if (r) {PetscValidHeaderSpecific(r,VEC_COOKIE);PetscCheckSameComm(mat,r);}
2347   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2348   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2349 
2350   PLogEventBegin(MAT_Scale,mat,0,0,0);
2351   ierr = (*mat->ops->diagonalscale)(mat,l,r);CHKERRQ(ierr);
2352   PLogEventEnd(MAT_Scale,mat,0,0,0);
2353   PetscFunctionReturn(0);
2354 }
2355 
2356 #undef __FUNC__
2357 #define __FUNC__ /*<a name=""></a>*/"MatScale"
2358 /*@
2359     MatScale - Scales all elements of a matrix by a given number.
2360 
2361     Collective on Mat
2362 
2363     Input Parameters:
2364 +   mat - the matrix to be scaled
2365 -   a  - the scaling value
2366 
2367     Output Parameter:
2368 .   mat - the scaled matrix
2369 
2370     Level: intermediate
2371 
2372 .keywords: matrix, scale
2373 
2374 .seealso: MatDiagonalScale()
2375 @*/
2376 int MatScale(Scalar *a,Mat mat)
2377 {
2378   int ierr;
2379 
2380   PetscFunctionBegin;
2381   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2382   PetscValidScalarPointer(a);
2383   if (!mat->ops->scale) SETERRQ(PETSC_ERR_SUP,0,"");
2384   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2385   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2386 
2387   PLogEventBegin(MAT_Scale,mat,0,0,0);
2388   ierr = (*mat->ops->scale)(a,mat);CHKERRQ(ierr);
2389   PLogEventEnd(MAT_Scale,mat,0,0,0);
2390   PetscFunctionReturn(0);
2391 }
2392 
2393 #undef __FUNC__
2394 #define __FUNC__ /*<a name=""></a>*/"MatNorm"
2395 /*@
2396    MatNorm - Calculates various norms of a matrix.
2397 
2398    Collective on Mat
2399 
2400    Input Parameters:
2401 +  mat - the matrix
2402 -  type - the type of norm, NORM_1, NORM_2, NORM_FROBENIUS, NORM_INFINITY
2403 
2404    Output Parameters:
2405 .  norm - the resulting norm
2406 
2407    Level: intermediate
2408 
2409 .keywords: matrix, norm, Frobenius
2410 @*/
2411 int MatNorm(Mat mat,NormType type,PetscReal *norm)
2412 {
2413   int ierr;
2414 
2415   PetscFunctionBegin;
2416   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2417   PetscValidScalarPointer(norm);
2418 
2419   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2420   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2421   if (!mat->ops->norm) SETERRQ(PETSC_ERR_SUP,0,"Not for this matrix type");
2422   ierr = (*mat->ops->norm)(mat,type,norm);CHKERRQ(ierr);
2423   PetscFunctionReturn(0);
2424 }
2425 
2426 /*
2427      This variable is used to prevent counting of MatAssemblyBegin() that
2428    are called from within a MatAssemblyEnd().
2429 */
2430 static int MatAssemblyEnd_InUse = 0;
2431 #undef __FUNC__
2432 #define __FUNC__ /*<a name=""></a>*/"MatAssemblyBegin"
2433 /*@
2434    MatAssemblyBegin - Begins assembling the matrix.  This routine should
2435    be called after completing all calls to MatSetValues().
2436 
2437    Collective on Mat
2438 
2439    Input Parameters:
2440 +  mat - the matrix
2441 -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
2442 
2443    Notes:
2444    MatSetValues() generally caches the values.  The matrix is ready to
2445    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
2446    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
2447    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
2448    using the matrix.
2449 
2450    Level: beginner
2451 
2452 .keywords: matrix, assembly, assemble, begin
2453 
2454 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled()
2455 @*/
2456 int MatAssemblyBegin(Mat mat,MatAssemblyType type)
2457 {
2458   int ierr;
2459 
2460   PetscFunctionBegin;
2461   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2462   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix.\n did you forget to call MatSetUnfactored()?");
2463   if (mat->assembled) {
2464     mat->was_assembled = PETSC_TRUE;
2465     mat->assembled     = PETSC_FALSE;
2466   }
2467   if (!MatAssemblyEnd_InUse) {
2468     PLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);
2469     if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);}
2470     PLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);
2471   } else {
2472     if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);}
2473   }
2474   PetscFunctionReturn(0);
2475 }
2476 
2477 #undef __FUNC__
2478 #define __FUNC__ /*<a name=""></a>*/"MatAssembed"
2479 /*@
2480    MatAssembled - Indicates if a matrix has been assembled and is ready for
2481      use; for example, in matrix-vector product.
2482 
2483    Collective on Mat
2484 
2485    Input Parameter:
2486 .  mat - the matrix
2487 
2488    Output Parameter:
2489 .  assembled - PETSC_TRUE or PETSC_FALSE
2490 
2491    Level: advanced
2492 
2493 .keywords: matrix, assembly, assemble, begin
2494 
2495 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin()
2496 @*/
2497 int MatAssembled(Mat mat,PetscTruth *assembled)
2498 {
2499   PetscFunctionBegin;
2500   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2501   *assembled = mat->assembled;
2502   PetscFunctionReturn(0);
2503 }
2504 
2505 #undef __FUNC__
2506 #define __FUNC__ /*<a name=""></a>*/"MatView_Private"
2507 /*
2508     Processes command line options to determine if/how a matrix
2509   is to be viewed. Called by MatAssemblyEnd() and MatLoad().
2510 */
2511 int MatView_Private(Mat mat)
2512 {
2513   int        ierr;
2514   PetscTruth flg;
2515 
2516   PetscFunctionBegin;
2517   ierr = OptionsHasName(mat->prefix,"-mat_view_info",&flg);CHKERRQ(ierr);
2518   if (flg) {
2519     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO,0);CHKERRQ(ierr);
2520     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2521     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2522   }
2523   ierr = OptionsHasName(mat->prefix,"-mat_view_info_detailed",&flg);CHKERRQ(ierr);
2524   if (flg) {
2525     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO_LONG,0);CHKERRQ(ierr);
2526     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2527     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2528   }
2529   ierr = OptionsHasName(mat->prefix,"-mat_view",&flg);CHKERRQ(ierr);
2530   if (flg) {
2531     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2532   }
2533   ierr = OptionsHasName(mat->prefix,"-mat_view_matlab",&flg);CHKERRQ(ierr);
2534   if (flg) {
2535     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_MATLAB,"M");CHKERRQ(ierr);
2536     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2537     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2538   }
2539   ierr = OptionsHasName(mat->prefix,"-mat_view_draw",&flg);CHKERRQ(ierr);
2540   if (flg) {
2541     ierr = OptionsHasName(mat->prefix,"-mat_view_contour",&flg);CHKERRQ(ierr);
2542     if (flg) {
2543       ViewerPushFormat(VIEWER_DRAW_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr);
2544     }
2545     ierr = MatView(mat,VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
2546     ierr = ViewerFlush(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
2547     if (flg) {
2548       ViewerPopFormat(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
2549     }
2550   }
2551   ierr = OptionsHasName(mat->prefix,"-mat_view_socket",&flg);CHKERRQ(ierr);
2552   if (flg) {
2553     ierr = MatView(mat,VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr);
2554     ierr = ViewerFlush(VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr);
2555   }
2556   PetscFunctionReturn(0);
2557 }
2558 
2559 #undef __FUNC__
2560 #define __FUNC__ /*<a name=""></a>*/"MatAssemblyEnd"
2561 /*@
2562    MatAssemblyEnd - Completes assembling the matrix.  This routine should
2563    be called after MatAssemblyBegin().
2564 
2565    Collective on Mat
2566 
2567    Input Parameters:
2568 +  mat - the matrix
2569 -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
2570 
2571    Options Database Keys:
2572 +  -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
2573 .  -mat_view_info_detailed - Prints more detailed info
2574 .  -mat_view - Prints matrix in ASCII format
2575 .  -mat_view_matlab - Prints matrix in Matlab format
2576 .  -mat_view_draw - Draws nonzero structure of matrix, using MatView() and DrawOpenX().
2577 .  -display <name> - Sets display name (default is host)
2578 -  -draw_pause <sec> - Sets number of seconds to pause after display
2579 
2580    Notes:
2581    MatSetValues() generally caches the values.  The matrix is ready to
2582    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
2583    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
2584    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
2585    using the matrix.
2586 
2587    Level: beginner
2588 
2589 .keywords: matrix, assembly, assemble, end
2590 
2591 .seealso: MatAssemblyBegin(), MatSetValues(), DrawOpenX(), MatView(), MatAssembled()
2592 @*/
2593 int MatAssemblyEnd(Mat mat,MatAssemblyType type)
2594 {
2595   int        ierr;
2596   static int inassm = 0;
2597 
2598   PetscFunctionBegin;
2599   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2600 
2601   inassm++;
2602   MatAssemblyEnd_InUse++;
2603   if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
2604     PLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);
2605     if (mat->ops->assemblyend) {
2606       ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr);
2607     }
2608     PLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);
2609   } else {
2610     if (mat->ops->assemblyend) {
2611       ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr);
2612     }
2613   }
2614 
2615   /* Flush assembly is not a true assembly */
2616   if (type != MAT_FLUSH_ASSEMBLY) {
2617     mat->assembled  = PETSC_TRUE; mat->num_ass++;
2618   }
2619   mat->insertmode = NOT_SET_VALUES;
2620   MatAssemblyEnd_InUse--;
2621 
2622   if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) {
2623     ierr = MatView_Private(mat);CHKERRQ(ierr);
2624   }
2625   inassm--;
2626   PetscFunctionReturn(0);
2627 }
2628 
2629 
2630 #undef __FUNC__
2631 #define __FUNC__ /*<a name=""></a>*/"MatCompress"
2632 /*@
2633    MatCompress - Tries to store the matrix in as little space as
2634    possible.  May fail if memory is already fully used, since it
2635    tries to allocate new space.
2636 
2637    Collective on Mat
2638 
2639    Input Parameters:
2640 .  mat - the matrix
2641 
2642    Level: advanced
2643 
2644 .keywords: matrix, compress
2645 @*/
2646 int MatCompress(Mat mat)
2647 {
2648   int ierr;
2649 
2650   PetscFunctionBegin;
2651   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2652   if (mat->ops->compress) {ierr = (*mat->ops->compress)(mat);CHKERRQ(ierr);}
2653   PetscFunctionReturn(0);
2654 }
2655 
2656 #undef __FUNC__
2657 #define __FUNC__ /*<a name=""></a>*/"MatSetOption"
2658 /*@
2659    MatSetOption - Sets a parameter option for a matrix. Some options
2660    may be specific to certain storage formats.  Some options
2661    determine how values will be inserted (or added). Sorted,
2662    row-oriented input will generally assemble the fastest. The default
2663    is row-oriented, nonsorted input.
2664 
2665    Collective on Mat
2666 
2667    Input Parameters:
2668 +  mat - the matrix
2669 -  option - the option, one of those listed below (and possibly others),
2670              e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR
2671 
2672    Options Describing Matrix Structure:
2673 +    MAT_SYMMETRIC - symmetric in terms of both structure and value
2674 -    MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure
2675 
2676    Options For Use with MatSetValues():
2677    Insert a logically dense subblock, which can be
2678 +    MAT_ROW_ORIENTED - row-oriented
2679 .    MAT_COLUMN_ORIENTED - column-oriented
2680 .    MAT_ROWS_SORTED - sorted by row
2681 .    MAT_ROWS_UNSORTED - not sorted by row
2682 .    MAT_COLUMNS_SORTED - sorted by column
2683 -    MAT_COLUMNS_UNSORTED - not sorted by column
2684 
2685    Not these options reflect the data you pass in with MatSetValues(); it has
2686    nothing to do with how the data is stored internally in the matrix
2687    data structure.
2688 
2689    When (re)assembling a matrix, we can restrict the input for
2690    efficiency/debugging purposes.  These options include
2691 +    MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be
2692         allowed if they generate a new nonzero
2693 .    MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed
2694 .    MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if
2695          they generate a nonzero in a new diagonal (for block diagonal format only)
2696 .    MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only)
2697 .    MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries
2698 .    MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry
2699 -    MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly
2700 
2701    Notes:
2702    Some options are relevant only for particular matrix types and
2703    are thus ignored by others.  Other options are not supported by
2704    certain matrix types and will generate an error message if set.
2705 
2706    If using a Fortran 77 module to compute a matrix, one may need to
2707    use the column-oriented option (or convert to the row-oriented
2708    format).
2709 
2710    MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion
2711    that would generate a new entry in the nonzero structure is instead
2712    ignored.  Thus, if memory has not alredy been allocated for this particular
2713    data, then the insertion is ignored. For dense matrices, in which
2714    the entire array is allocated, no entries are ever ignored.
2715 
2716    MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion
2717    that would generate a new entry in the nonzero structure instead produces
2718    an error. (Currently supported for AIJ and BAIJ formats only.)
2719    This is a useful flag when using SAME_NONZERO_PATTERN in calling
2720    SLESSetOperators() to ensure that the nonzero pattern truely does
2721    remain unchanged.
2722 
2723    MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion
2724    that would generate a new entry that has not been preallocated will
2725    instead produce an error. (Currently supported for AIJ and BAIJ formats
2726    only.) This is a useful flag when debugging matrix memory preallocation.
2727 
2728    MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for
2729    other processors should be dropped, rather than stashed.
2730    This is useful if you know that the "owning" processor is also
2731    always generating the correct matrix entries, so that PETSc need
2732    not transfer duplicate entries generated on another processor.
2733 
2734    MAT_USE_HASH_TABLE indicates that a hash table be used to improve the
2735    searches during matrix assembly. When this flag is set, the hash table
2736    is created during the first Matrix Assembly. This hash table is
2737    used the next time through, during MatSetVaules()/MatSetVaulesBlocked()
2738    to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag
2739    should be used with MAT_USE_HASH_TABLE flag. This option is currently
2740    supported by MATMPIBAIJ format only.
2741 
2742    MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries
2743    are kept in the nonzero structure
2744 
2745    MAT_IGNORE_ZERO_ENTRIES - when using ADD_VALUES for AIJ matrices this will stop
2746    zero values from creating a zero location in the matrix
2747 
2748    Level: intermediate
2749 
2750 .keywords: matrix, option, row-oriented, column-oriented, sorted, nonzero
2751 @*/
2752 int MatSetOption(Mat mat,MatOption op)
2753 {
2754   int ierr;
2755 
2756   PetscFunctionBegin;
2757   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2758   if (mat->ops->setoption) {ierr = (*mat->ops->setoption)(mat,op);CHKERRQ(ierr);}
2759   PetscFunctionReturn(0);
2760 }
2761 
2762 #undef __FUNC__
2763 #define __FUNC__ /*<a name=""></a>*/"MatZeroEntries"
2764 /*@
2765    MatZeroEntries - Zeros all entries of a matrix.  For sparse matrices
2766    this routine retains the old nonzero structure.
2767 
2768    Collective on Mat
2769 
2770    Input Parameters:
2771 .  mat - the matrix
2772 
2773    Level: intermediate
2774 
2775 .keywords: matrix, zero, entries
2776 
2777 .seealso: MatZeroRows()
2778 @*/
2779 int MatZeroEntries(Mat mat)
2780 {
2781   int ierr;
2782 
2783   PetscFunctionBegin;
2784   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2785   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2786   if (!mat->ops->zeroentries) SETERRQ(PETSC_ERR_SUP,0,"");
2787 
2788   PLogEventBegin(MAT_ZeroEntries,mat,0,0,0);
2789   ierr = (*mat->ops->zeroentries)(mat);CHKERRQ(ierr);
2790   PLogEventEnd(MAT_ZeroEntries,mat,0,0,0);
2791   PetscFunctionReturn(0);
2792 }
2793 
2794 #undef __FUNC__
2795 #define __FUNC__ /*<a name=""></a>*/"MatZeroRows"
2796 /*@C
2797    MatZeroRows - Zeros all entries (except possibly the main diagonal)
2798    of a set of rows of a matrix.
2799 
2800    Collective on Mat
2801 
2802    Input Parameters:
2803 +  mat - the matrix
2804 .  is - index set of rows to remove
2805 -  diag - pointer to value put in all diagonals of eliminated rows.
2806           Note that diag is not a pointer to an array, but merely a
2807           pointer to a single value.
2808 
2809    Notes:
2810    For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
2811    but does not release memory.  For the dense and block diagonal
2812    formats this does not alter the nonzero structure.
2813 
2814    If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
2815    of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
2816    merely zeroed.
2817 
2818    The user can set a value in the diagonal entry (or for the AIJ and
2819    row formats can optionally remove the main diagonal entry from the
2820    nonzero structure as well, by passing a null pointer (PETSC_NULL
2821    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
2822 
2823    For the parallel case, all processes that share the matrix (i.e.,
2824    those in the communicator used for matrix creation) MUST call this
2825    routine, regardless of whether any rows being zeroed are owned by
2826    them.
2827 
2828 
2829    Level: intermediate
2830 
2831 .keywords: matrix, zero, rows, boundary conditions
2832 
2833 .seealso: MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
2834 @*/
2835 int MatZeroRows(Mat mat,IS is,Scalar *diag)
2836 {
2837   int ierr;
2838 
2839   PetscFunctionBegin;
2840   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2841   PetscValidHeaderSpecific(is,IS_COOKIE);
2842   if (diag) PetscValidScalarPointer(diag);
2843   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2844   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2845   if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,0,"");
2846 
2847   ierr = (*mat->ops->zerorows)(mat,is,diag);CHKERRQ(ierr);
2848   ierr = MatView_Private(mat);CHKERRQ(ierr);
2849   PetscFunctionReturn(0);
2850 }
2851 
2852 #undef __FUNC__
2853 #define __FUNC__ /*<a name=""></a>*/"MatZeroRowsLocal"
2854 /*@C
2855    MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
2856    of a set of rows of a matrix; using local numbering of rows.
2857 
2858    Collective on Mat
2859 
2860    Input Parameters:
2861 +  mat - the matrix
2862 .  is - index set of rows to remove
2863 -  diag - pointer to value put in all diagonals of eliminated rows.
2864           Note that diag is not a pointer to an array, but merely a
2865           pointer to a single value.
2866 
2867    Notes:
2868    For the AIJ matrix formats this removes the old nonzero structure,
2869    but does not release memory.  For the dense and block diagonal
2870    formats this does not alter the nonzero structure.
2871 
2872    The user can set a value in the diagonal entry (or for the AIJ and
2873    row formats can optionally remove the main diagonal entry from the
2874    nonzero structure as well, by passing a null pointer (PETSC_NULL
2875    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
2876 
2877    Level: intermediate
2878 
2879 .keywords: matrix, zero, rows, boundary conditions
2880 
2881 .seealso: MatZeroEntries(), MatZeroRows()
2882 @*/
2883 int MatZeroRowsLocal(Mat mat,IS is,Scalar *diag)
2884 {
2885   int ierr;
2886   IS  newis;
2887 
2888   PetscFunctionBegin;
2889   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2890   PetscValidHeaderSpecific(is,IS_COOKIE);
2891   if (diag) PetscValidScalarPointer(diag);
2892   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2893   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2894   if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,0,"");
2895 
2896   ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);CHKERRQ(ierr);
2897   ierr = (*mat->ops->zerorows)(mat,newis,diag);CHKERRQ(ierr);
2898   ierr = ISDestroy(newis);CHKERRQ(ierr);
2899   PetscFunctionReturn(0);
2900 }
2901 
2902 #undef __FUNC__
2903 #define __FUNC__ /*<a name=""></a>*/"MatGetSize"
2904 /*@
2905    MatGetSize - Returns the numbers of rows and columns in a matrix.
2906 
2907    Not Collective
2908 
2909    Input Parameter:
2910 .  mat - the matrix
2911 
2912    Output Parameters:
2913 +  m - the number of global rows
2914 -  n - the number of global columns
2915 
2916    Level: beginner
2917 
2918 .keywords: matrix, dimension, size, rows, columns, global, get
2919 
2920 .seealso: MatGetLocalSize()
2921 @*/
2922 int MatGetSize(Mat mat,int *m,int* n)
2923 {
2924   int ierr;
2925 
2926   PetscFunctionBegin;
2927   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2928   ierr = (*mat->ops->getsize)(mat,m,n);CHKERRQ(ierr);
2929   PetscFunctionReturn(0);
2930 }
2931 
2932 #undef __FUNC__
2933 #define __FUNC__ /*<a name=""></a>*/"MatGetLocalSize"
2934 /*@
2935    MatGetLocalSize - Returns the number of rows and columns in a matrix
2936    stored locally.  This information may be implementation dependent, so
2937    use with care.
2938 
2939    Not Collective
2940 
2941    Input Parameters:
2942 .  mat - the matrix
2943 
2944    Output Parameters:
2945 +  m - the number of local rows
2946 -  n - the number of local columns
2947 
2948    Level: beginner
2949 
2950 .keywords: matrix, dimension, size, local, rows, columns, get
2951 
2952 .seealso: MatGetSize()
2953 @*/
2954 int MatGetLocalSize(Mat mat,int *m,int* n)
2955 {
2956   int ierr;
2957 
2958   PetscFunctionBegin;
2959   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2960   ierr = (*mat->ops->getlocalsize)(mat,m,n);CHKERRQ(ierr);
2961   PetscFunctionReturn(0);
2962 }
2963 
2964 #undef __FUNC__
2965 #define __FUNC__ /*<a name=""></a>*/"MatGetOwnershipRange"
2966 /*@
2967    MatGetOwnershipRange - Returns the range of matrix rows owned by
2968    this processor, assuming that the matrix is laid out with the first
2969    n1 rows on the first processor, the next n2 rows on the second, etc.
2970    For certain parallel layouts this range may not be well defined.
2971 
2972    Not Collective
2973 
2974    Input Parameters:
2975 .  mat - the matrix
2976 
2977    Output Parameters:
2978 +  m - the global index of the first local row
2979 -  n - one more than the global index of the last local row
2980 
2981    Level: beginner
2982 
2983 .keywords: matrix, get, range, ownership
2984 @*/
2985 int MatGetOwnershipRange(Mat mat,int *m,int* n)
2986 {
2987   int ierr;
2988 
2989   PetscFunctionBegin;
2990   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2991   if (m) PetscValidIntPointer(m);
2992   if (n) PetscValidIntPointer(n);
2993   if (!mat->ops->getownershiprange) SETERRQ(PETSC_ERR_SUP,0,"");
2994   ierr = (*mat->ops->getownershiprange)(mat,m,n);CHKERRQ(ierr);
2995   PetscFunctionReturn(0);
2996 }
2997 
2998 #undef __FUNC__
2999 #define __FUNC__ /*<a name=""></a>*/"MatILUFactorSymbolic"
3000 /*@
3001    MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
3002    Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric()
3003    to complete the factorization.
3004 
3005    Collective on Mat
3006 
3007    Input Parameters:
3008 +  mat - the matrix
3009 .  row - row permutation
3010 .  column - column permutation
3011 -  info - structure containing
3012 $      levels - number of levels of fill.
3013 $      expected fill - as ratio of original fill.
3014 $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
3015                 missing diagonal entries)
3016 
3017    Output Parameters:
3018 .  fact - new matrix that has been symbolically factored
3019 
3020    Notes:
3021    See the users manual for additional information about
3022    choosing the fill factor for better efficiency.
3023 
3024    Most users should employ the simplified SLES interface for linear solvers
3025    instead of working directly with matrix algebra routines such as this.
3026    See, e.g., SLESCreate().
3027 
3028    Level: developer
3029 
3030 .keywords: matrix, factor, incomplete, ILU, symbolic, fill
3031 
3032 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
3033           MatGetOrdering()
3034 
3035 @*/
3036 int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatILUInfo *info,Mat *fact)
3037 {
3038   int ierr;
3039 
3040   PetscFunctionBegin;
3041   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3042   PetscValidPointer(fact);
3043   if (info && info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Levels of fill negative %d",info->levels);
3044   if (info && info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",info->fill);
3045   if (!mat->ops->ilufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Only MatCreateMPIRowbs() matrices support parallel ILU");
3046   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3047   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3048 
3049   PLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);
3050   ierr = (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr);
3051   PLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);
3052   PetscFunctionReturn(0);
3053 }
3054 
3055 #undef __FUNC__
3056 #define __FUNC__ /*<a name=""></a>*/"MatIncompleteCholeskyFactorSymbolic"
3057 /*@
3058    MatIncompleteCholeskyFactorSymbolic - Performs symbolic incomplete
3059    Cholesky factorization for a symmetric matrix.  Use
3060    MatCholeskyFactorNumeric() to complete the factorization.
3061 
3062    Collective on Mat
3063 
3064    Input Parameters:
3065 +  mat - the matrix
3066 .  perm - row and column permutation
3067 .  fill - levels of fill
3068 -  f - expected fill as ratio of original fill
3069 
3070    Output Parameter:
3071 .  fact - the factored matrix
3072 
3073    Notes:
3074    Currently only no-fill factorization is supported.
3075 
3076    Most users should employ the simplified SLES interface for linear solvers
3077    instead of working directly with matrix algebra routines such as this.
3078    See, e.g., SLESCreate().
3079 
3080    Level: developer
3081 
3082 .keywords: matrix, factor, incomplete, ICC, Cholesky, symbolic, fill
3083 
3084 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor()
3085 @*/
3086 int MatIncompleteCholeskyFactorSymbolic(Mat mat,IS perm,PetscReal f,int fill,Mat *fact)
3087 {
3088   int ierr;
3089 
3090   PetscFunctionBegin;
3091   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3092   PetscValidPointer(fact);
3093   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3094   if (fill < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Fill negative %d",fill);
3095   if (f < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",f);
3096   if (!mat->ops->incompletecholeskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support ICC in parallel");
3097   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3098 
3099   PLogEventBegin(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);
3100   ierr = (*mat->ops->incompletecholeskyfactorsymbolic)(mat,perm,f,fill,fact);CHKERRQ(ierr);
3101   PLogEventEnd(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);
3102   PetscFunctionReturn(0);
3103 }
3104 
3105 #undef __FUNC__
3106 #define __FUNC__ /*<a name=""></a>*/"MatGetArray"
3107 /*@C
3108    MatGetArray - Returns a pointer to the element values in the matrix.
3109    The result of this routine is dependent on the underlying matrix data
3110    structure, and may not even work for certain matrix types.  You MUST
3111    call MatRestoreArray() when you no longer need to access the array.
3112 
3113    Not Collective
3114 
3115    Input Parameter:
3116 .  mat - the matrix
3117 
3118    Output Parameter:
3119 .  v - the location of the values
3120 
3121    Currently returns an array only for the dense formats, giving access to
3122    the local portion of the matrix in the usual Fortran column-oriented format.
3123 
3124    Fortran Note:
3125    This routine is used differently from Fortran, e.g.,
3126 .vb
3127         Mat         mat
3128         Scalar      mat_array(1)
3129         PetscOffset i_mat
3130         int         ierr
3131         call MatGetArray(mat,mat_array,i_mat,ierr)
3132 
3133   C  Access first local entry in matrix; note that array is
3134   C  treated as one dimensional
3135         value = mat_array(i_mat + 1)
3136 
3137         [... other code ...]
3138         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3139 .ve
3140 
3141    See the Fortran chapter of the users manual and
3142    petsc/src/mat/examples/tests for details.
3143 
3144    Level: advanced
3145 
3146 .keywords: matrix, array, elements, values
3147 
3148 .seealso: MatRestoreArray(), MatGetArrayF90()
3149 @*/
3150 int MatGetArray(Mat mat,Scalar **v)
3151 {
3152   int ierr;
3153 
3154   PetscFunctionBegin;
3155   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3156   PetscValidPointer(v);
3157   if (!mat->ops->getarray) SETERRQ(PETSC_ERR_SUP,0,"");
3158   ierr = (*mat->ops->getarray)(mat,v);CHKERRQ(ierr);
3159   PetscFunctionReturn(0);
3160 }
3161 
3162 #undef __FUNC__
3163 #define __FUNC__ /*<a name=""></a>*/"MatRestoreArray"
3164 /*@C
3165    MatRestoreArray - Restores the matrix after MatGetArray() has been called.
3166 
3167    Not Collective
3168 
3169    Input Parameter:
3170 +  mat - the matrix
3171 -  v - the location of the values
3172 
3173    Fortran Note:
3174    This routine is used differently from Fortran, e.g.,
3175 .vb
3176         Mat         mat
3177         Scalar      mat_array(1)
3178         PetscOffset i_mat
3179         int         ierr
3180         call MatGetArray(mat,mat_array,i_mat,ierr)
3181 
3182   C  Access first local entry in matrix; note that array is
3183   C  treated as one dimensional
3184         value = mat_array(i_mat + 1)
3185 
3186         [... other code ...]
3187         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3188 .ve
3189 
3190    See the Fortran chapter of the users manual and
3191    petsc/src/mat/examples/tests for details
3192 
3193    Level: advanced
3194 
3195 .keywords: matrix, array, elements, values, restore
3196 
3197 .seealso: MatGetArray(), MatRestoreArrayF90()
3198 @*/
3199 int MatRestoreArray(Mat mat,Scalar **v)
3200 {
3201   int ierr;
3202 
3203   PetscFunctionBegin;
3204   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3205   PetscValidPointer(v);
3206   if (!mat->ops->restorearray) SETERRQ(PETSC_ERR_SUP,0,"");
3207   ierr = (*mat->ops->restorearray)(mat,v);CHKERRQ(ierr);
3208   PetscFunctionReturn(0);
3209 }
3210 
3211 #undef __FUNC__
3212 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrices"
3213 /*@C
3214    MatGetSubMatrices - Extracts several submatrices from a matrix. If submat
3215    points to an array of valid matrices, they may be reused to store the new
3216    submatrices.
3217 
3218    Collective on Mat
3219 
3220    Input Parameters:
3221 +  mat - the matrix
3222 .  n   - the number of submatrixes to be extracted (on this processor, may be zero)
3223 .  irow, icol - index sets of rows and columns to extract
3224 -  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3225 
3226    Output Parameter:
3227 .  submat - the array of submatrices
3228 
3229    Notes:
3230    MatGetSubMatrices() can extract only sequential submatrices
3231    (from both sequential and parallel matrices). Use MatGetSubMatrix()
3232    to extract a parallel submatrix.
3233 
3234    When extracting submatrices from a parallel matrix, each processor can
3235    form a different submatrix by setting the rows and columns of its
3236    individual index sets according to the local submatrix desired.
3237 
3238    When finished using the submatrices, the user should destroy
3239    them with MatDestroySubMatrices().
3240 
3241    MAT_REUSE_MATRIX can only be used when the nonzero structure of the
3242    original matrix has not changed from that last call to MatGetSubMatrices().
3243 
3244    Fortran Note:
3245    The Fortran interface is slightly different from that given below; it
3246    requires one to pass in  as submat a Mat (integer) array of size at least m.
3247 
3248    Level: advanced
3249 
3250 .keywords: matrix, get, submatrix, submatrices
3251 
3252 .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal()
3253 @*/
3254 int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatReuse scall,Mat **submat)
3255 {
3256   int ierr;
3257 
3258   PetscFunctionBegin;
3259   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3260   if (!mat->ops->getsubmatrices) SETERRQ(PETSC_ERR_SUP,0,"");
3261   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3262 
3263   PLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);
3264   ierr = (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);CHKERRQ(ierr);
3265   PLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);
3266 
3267   PetscFunctionReturn(0);
3268 }
3269 
3270 #undef __FUNC__
3271 #define __FUNC__ /*<a name=""></a>*/"MatDestroyMatrices"
3272 /*@C
3273    MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices().
3274 
3275    Collective on Mat
3276 
3277    Input Parameters:
3278 +  n - the number of local matrices
3279 -  mat - the matrices
3280 
3281    Level: advanced
3282 
3283 .keywords: matrix, destroy, submatrix, submatrices
3284 
3285 .seealso: MatGetSubMatrices()
3286 @*/
3287 int MatDestroyMatrices(int n,Mat **mat)
3288 {
3289   int ierr,i;
3290 
3291   PetscFunctionBegin;
3292   if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,1,"Trying to destroy negative number of matrices %d",n);
3293   PetscValidPointer(mat);
3294   for (i=0; i<n; i++) {
3295     ierr = MatDestroy((*mat)[i]);CHKERRQ(ierr);
3296   }
3297   if (n) {ierr = PetscFree(*mat);CHKERRQ(ierr);}
3298   PetscFunctionReturn(0);
3299 }
3300 
3301 #undef __FUNC__
3302 #define __FUNC__ /*<a name=""></a>*/"MatIncreaseOverlap"
3303 /*@
3304    MatIncreaseOverlap - Given a set of submatrices indicated by index sets,
3305    replaces the index sets by larger ones that represent submatrices with
3306    additional overlap.
3307 
3308    Collective on Mat
3309 
3310    Input Parameters:
3311 +  mat - the matrix
3312 .  n   - the number of index sets
3313 .  is  - the array of pointers to index sets
3314 -  ov  - the additional overlap requested
3315 
3316    Level: developer
3317 
3318 .keywords: matrix, overlap, Schwarz
3319 
3320 .seealso: MatGetSubMatrices()
3321 @*/
3322 int MatIncreaseOverlap(Mat mat,int n,IS *is,int ov)
3323 {
3324   int ierr;
3325 
3326   PetscFunctionBegin;
3327   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3328   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3329   if (mat->factor)     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3330 
3331   if (!ov) PetscFunctionReturn(0);
3332   if (!mat->ops->increaseoverlap) SETERRQ(PETSC_ERR_SUP,0,"");
3333   PLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);
3334   ierr = (*mat->ops->increaseoverlap)(mat,n,is,ov);CHKERRQ(ierr);
3335   PLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);
3336   PetscFunctionReturn(0);
3337 }
3338 
3339 #undef __FUNC__
3340 #define __FUNC__ /*<a name=""></a>*/"MatPrintHelp"
3341 /*@
3342    MatPrintHelp - Prints all the options for the matrix.
3343 
3344    Collective on Mat
3345 
3346    Input Parameter:
3347 .  mat - the matrix
3348 
3349    Options Database Keys:
3350 +  -help - Prints matrix options
3351 -  -h - Prints matrix options
3352 
3353    Level: developer
3354 
3355 .keywords: mat, help
3356 
3357 .seealso: MatCreate(), MatCreateXXX()
3358 @*/
3359 int MatPrintHelp(Mat mat)
3360 {
3361   static PetscTruth called = PETSC_FALSE;
3362   int               ierr;
3363   MPI_Comm          comm;
3364 
3365   PetscFunctionBegin;
3366   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3367 
3368   comm = mat->comm;
3369   if (!called) {
3370     ierr = (*PetscHelpPrintf)(comm,"General matrix options:\n");CHKERRQ(ierr);
3371     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_info: view basic matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr);
3372     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_info_detailed: view detailed matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr);
3373     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_draw: draw nonzero matrix structure during MatAssemblyEnd()\n");CHKERRQ(ierr);
3374     ierr = (*PetscHelpPrintf)(comm,"      -draw_pause <sec>: set seconds of display pause\n");CHKERRQ(ierr);
3375     ierr = (*PetscHelpPrintf)(comm,"      -display <name>: set alternate display\n");CHKERRQ(ierr);
3376     called = PETSC_TRUE;
3377   }
3378   if (mat->ops->printhelp) {
3379     ierr = (*mat->ops->printhelp)(mat);CHKERRQ(ierr);
3380   }
3381   PetscFunctionReturn(0);
3382 }
3383 
3384 #undef __FUNC__
3385 #define __FUNC__ /*<a name=""></a>*/"MatGetBlockSize"
3386 /*@
3387    MatGetBlockSize - Returns the matrix block size; useful especially for the
3388    block row and block diagonal formats.
3389 
3390    Not Collective
3391 
3392    Input Parameter:
3393 .  mat - the matrix
3394 
3395    Output Parameter:
3396 .  bs - block size
3397 
3398    Notes:
3399    Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG.
3400    Block row formats are MATSEQBAIJ, MATMPIBAIJ
3401 
3402    Level: intermediate
3403 
3404 .keywords: matrix, get, block, size
3405 
3406 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag()
3407 @*/
3408 int MatGetBlockSize(Mat mat,int *bs)
3409 {
3410   int ierr;
3411 
3412   PetscFunctionBegin;
3413   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3414   PetscValidIntPointer(bs);
3415   if (!mat->ops->getblocksize) SETERRQ(PETSC_ERR_SUP,0,"");
3416   ierr = (*mat->ops->getblocksize)(mat,bs);CHKERRQ(ierr);
3417   PetscFunctionReturn(0);
3418 }
3419 
3420 #undef __FUNC__
3421 #define __FUNC__ /*<a name=""></a>*/"MatGetRowIJ"
3422 /*@C
3423     MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices.
3424 
3425    Collective on Mat
3426 
3427     Input Parameters:
3428 +   mat - the matrix
3429 .   shift -  0 or 1 indicating we want the indices starting at 0 or 1
3430 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3431                 symmetrized
3432 
3433     Output Parameters:
3434 +   n - number of rows in the (possibly compressed) matrix
3435 .   ia - the row pointers
3436 .   ja - the column indices
3437 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3438 
3439     Level: developer
3440 
3441 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3442 @*/
3443 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3444 {
3445   int ierr;
3446 
3447   PetscFunctionBegin;
3448   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3449   if (ia) PetscValidIntPointer(ia);
3450   if (ja) PetscValidIntPointer(ja);
3451   PetscValidIntPointer(done);
3452   if (!mat->ops->getrowij) *done = PETSC_FALSE;
3453   else {
3454     *done = PETSC_TRUE;
3455     ierr  = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3456   }
3457   PetscFunctionReturn(0);
3458 }
3459 
3460 #undef __FUNC__
3461 #define __FUNC__ /*<a name=""></a>*/"MatGetColumnIJ"
3462 /*@C
3463     MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices.
3464 
3465     Collective on Mat
3466 
3467     Input Parameters:
3468 +   mat - the matrix
3469 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3470 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3471                 symmetrized
3472 
3473     Output Parameters:
3474 +   n - number of columns in the (possibly compressed) matrix
3475 .   ia - the column pointers
3476 .   ja - the row indices
3477 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3478 
3479     Level: developer
3480 
3481 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3482 @*/
3483 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3484 {
3485   int ierr;
3486 
3487   PetscFunctionBegin;
3488   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3489   if (ia) PetscValidIntPointer(ia);
3490   if (ja) PetscValidIntPointer(ja);
3491   PetscValidIntPointer(done);
3492 
3493   if (!mat->ops->getcolumnij) *done = PETSC_FALSE;
3494   else {
3495     *done = PETSC_TRUE;
3496     ierr  = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3497   }
3498   PetscFunctionReturn(0);
3499 }
3500 
3501 #undef __FUNC__
3502 #define __FUNC__ /*<a name=""></a>*/"MatRestoreRowIJ"
3503 /*@C
3504     MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with
3505     MatGetRowIJ().
3506 
3507     Collective on Mat
3508 
3509     Input Parameters:
3510 +   mat - the matrix
3511 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3512 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3513                 symmetrized
3514 
3515     Output Parameters:
3516 +   n - size of (possibly compressed) matrix
3517 .   ia - the row pointers
3518 .   ja - the column indices
3519 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3520 
3521     Level: developer
3522 
3523 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3524 @*/
3525 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3526 {
3527   int ierr;
3528 
3529   PetscFunctionBegin;
3530   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3531   if (ia) PetscValidIntPointer(ia);
3532   if (ja) PetscValidIntPointer(ja);
3533   PetscValidIntPointer(done);
3534 
3535   if (!mat->ops->restorerowij) *done = PETSC_FALSE;
3536   else {
3537     *done = PETSC_TRUE;
3538     ierr  = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3539   }
3540   PetscFunctionReturn(0);
3541 }
3542 
3543 #undef __FUNC__
3544 #define __FUNC__ /*<a name=""></a>*/"MatRestoreColumnIJ"
3545 /*@C
3546     MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with
3547     MatGetColumnIJ().
3548 
3549     Collective on Mat
3550 
3551     Input Parameters:
3552 +   mat - the matrix
3553 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3554 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3555                 symmetrized
3556 
3557     Output Parameters:
3558 +   n - size of (possibly compressed) matrix
3559 .   ia - the column pointers
3560 .   ja - the row indices
3561 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3562 
3563     Level: developer
3564 
3565 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3566 @*/
3567 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3568 {
3569   int ierr;
3570 
3571   PetscFunctionBegin;
3572   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3573   if (ia) PetscValidIntPointer(ia);
3574   if (ja) PetscValidIntPointer(ja);
3575   PetscValidIntPointer(done);
3576 
3577   if (!mat->ops->restorecolumnij) *done = PETSC_FALSE;
3578   else {
3579     *done = PETSC_TRUE;
3580     ierr  = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3581   }
3582   PetscFunctionReturn(0);
3583 }
3584 
3585 #undef __FUNC__
3586 #define __FUNC__ /*<a name=""></a>*/"MatColoringPatch"
3587 /*@C
3588     MatColoringPatch -Used inside matrix coloring routines that
3589     use MatGetRowIJ() and/or MatGetColumnIJ().
3590 
3591     Collective on Mat
3592 
3593     Input Parameters:
3594 +   mat - the matrix
3595 .   n   - number of colors
3596 -   colorarray - array indicating color for each column
3597 
3598     Output Parameters:
3599 .   iscoloring - coloring generated using colorarray information
3600 
3601     Level: developer
3602 
3603 .seealso: MatGetRowIJ(), MatGetColumnIJ()
3604 
3605 .keywords: mat, coloring, patch
3606 @*/
3607 int MatColoringPatch(Mat mat,int n,int *colorarray,ISColoring *iscoloring)
3608 {
3609   int ierr;
3610 
3611   PetscFunctionBegin;
3612   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3613   PetscValidIntPointer(colorarray);
3614 
3615   if (!mat->ops->coloringpatch) {SETERRQ(PETSC_ERR_SUP,0,"");}
3616   else {
3617     ierr  = (*mat->ops->coloringpatch)(mat,n,colorarray,iscoloring);CHKERRQ(ierr);
3618   }
3619   PetscFunctionReturn(0);
3620 }
3621 
3622 
3623 #undef __FUNC__
3624 #define __FUNC__ /*<a name=""></a>*/"MatSetUnfactored"
3625 /*@
3626    MatSetUnfactored - Resets a factored matrix to be treated as unfactored.
3627 
3628    Collective on Mat
3629 
3630    Input Parameter:
3631 .  mat - the factored matrix to be reset
3632 
3633    Notes:
3634    This routine should be used only with factored matrices formed by in-place
3635    factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE
3636    format).  This option can save memory, for example, when solving nonlinear
3637    systems with a matrix-free Newton-Krylov method and a matrix-based, in-place
3638    ILU(0) preconditioner.
3639 
3640    Note that one can specify in-place ILU(0) factorization by calling
3641 .vb
3642      PCType(pc,PCILU);
3643      PCILUSeUseInPlace(pc);
3644 .ve
3645    or by using the options -pc_type ilu -pc_ilu_in_place
3646 
3647    In-place factorization ILU(0) can also be used as a local
3648    solver for the blocks within the block Jacobi or additive Schwarz
3649    methods (runtime option: -sub_pc_ilu_in_place).  See the discussion
3650    of these preconditioners in the users manual for details on setting
3651    local solver options.
3652 
3653    Most users should employ the simplified SLES interface for linear solvers
3654    instead of working directly with matrix algebra routines such as this.
3655    See, e.g., SLESCreate().
3656 
3657    Level: developer
3658 
3659 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace()
3660 
3661 .keywords: matrix-free, in-place ILU, in-place LU
3662 @*/
3663 int MatSetUnfactored(Mat mat)
3664 {
3665   int ierr;
3666 
3667   PetscFunctionBegin;
3668   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3669   mat->factor = 0;
3670   if (!mat->ops->setunfactored) PetscFunctionReturn(0);
3671   ierr = (*mat->ops->setunfactored)(mat);CHKERRQ(ierr);
3672   PetscFunctionReturn(0);
3673 }
3674 
3675 #undef __FUNC__
3676 #define __FUNC__ /*<a name=""></a>*/"MatGetType"
3677 /*@C
3678    MatGetType - Gets the matrix type and name (as a string) from the matrix.
3679 
3680    Not Collective
3681 
3682    Input Parameter:
3683 .  mat - the matrix
3684 
3685    Output Parameter:
3686 +  type - the matrix type (or use PETSC_NULL)
3687 -  name - name of matrix type (or use PETSC_NULL)
3688 
3689    Level: intermediate
3690 
3691 .keywords: matrix, get, type, name
3692 @*/
3693 int MatGetType(Mat mat,MatType *type,char **name)
3694 {
3695   int  itype = (int)mat->type;
3696   char *matname[10];
3697 
3698   PetscFunctionBegin;
3699   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3700 
3701   if (type) *type = (MatType) mat->type;
3702   if (name) {
3703     /* Note:  Be sure that this list corresponds to the enum in petscmat.h */
3704     matname[0] = "MATSEQDENSE";
3705     matname[1] = "MATSEQAIJ";
3706     matname[2] = "MATMPIAIJ";
3707     matname[3] = "MATSHELL";
3708     matname[4] = "MATMPIROWBS";
3709     matname[5] = "MATSEQBDIAG";
3710     matname[6] = "MATMPIBDIAG";
3711     matname[7] = "MATMPIDENSE";
3712     matname[8] = "MATSEQBAIJ";
3713     matname[9] = "MATMPIBAIJ";
3714 
3715     if (itype < 0 || itype > 9) *name = "Unknown matrix type";
3716     else                        *name = matname[itype];
3717   }
3718   PetscFunctionReturn(0);
3719 }
3720 
3721 /*MC
3722     MatGetArrayF90 - Accesses a matrix array from Fortran90.
3723 
3724     Synopsis:
3725     MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3726 
3727     Not collective
3728 
3729     Input Parameter:
3730 .   x - matrix
3731 
3732     Output Parameters:
3733 +   xx_v - the Fortran90 pointer to the array
3734 -   ierr - error code
3735 
3736     Example of Usage:
3737 .vb
3738       Scalar, pointer xx_v(:)
3739       ....
3740       call MatGetArrayF90(x,xx_v,ierr)
3741       a = xx_v(3)
3742       call MatRestoreArrayF90(x,xx_v,ierr)
3743 .ve
3744 
3745     Notes:
3746     Not yet supported for all F90 compilers
3747 
3748     Level: advanced
3749 
3750 .seealso:  MatRestoreArrayF90(), MatGetArray(), MatRestoreArray()
3751 
3752 .keywords:  matrix, array, f90
3753 M*/
3754 
3755 /*MC
3756     MatRestoreArrayF90 - Restores a matrix array that has been
3757     accessed with MatGetArrayF90().
3758 
3759     Synopsis:
3760     MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3761 
3762     Not collective
3763 
3764     Input Parameters:
3765 +   x - matrix
3766 -   xx_v - the Fortran90 pointer to the array
3767 
3768     Output Parameter:
3769 .   ierr - error code
3770 
3771     Example of Usage:
3772 .vb
3773        Scalar, pointer xx_v(:)
3774        ....
3775        call MatGetArrayF90(x,xx_v,ierr)
3776        a = xx_v(3)
3777        call MatRestoreArrayF90(x,xx_v,ierr)
3778 .ve
3779 
3780     Notes:
3781     Not yet supported for all F90 compilers
3782 
3783     Level: advanced
3784 
3785 .seealso:  MatGetArrayF90(), MatGetArray(), MatRestoreArray()
3786 
3787 .keywords:  matrix, array, f90
3788 M*/
3789 
3790 
3791 #undef __FUNC__
3792 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrix"
3793 /*@
3794     MatGetSubMatrix - Gets a single submatrix on the same number of processors
3795                       as the original matrix.
3796 
3797     Collective on Mat
3798 
3799     Input Parameters:
3800 +   mat - the original matrix
3801 .   isrow - rows this processor should obtain
3802 .   iscol - columns for all processors you wish to keep
3803 .   csize - number of columns "local" to this processor (does nothing for sequential
3804             matrices). This should match the result from VecGetLocalSize(x,...) if you
3805             plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE
3806 -   cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3807 
3808     Output Parameter:
3809 .   newmat - the new submatrix, of the same type as the old
3810 
3811     Level: advanced
3812 
3813 .keywords: matrix, get, submatrix, submatrices
3814 
3815 .seealso: MatGetSubMatrices()
3816 @*/
3817 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat)
3818 {
3819   int     ierr, size;
3820   Mat     *local;
3821 
3822   PetscFunctionBegin;
3823   ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr);
3824 
3825   /* if original matrix is on just one processor then use submatrix generated */
3826   if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) {
3827     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr);
3828     PetscFunctionReturn(0);
3829   } else if (!mat->ops->getsubmatrix && size == 1) {
3830     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr);
3831     *newmat = *local;
3832     ierr = PetscFree(local);CHKERRQ(ierr);
3833     PetscFunctionReturn(0);
3834   }
3835 
3836   if (!mat->ops->getsubmatrix) SETERRQ(PETSC_ERR_SUP,0,"Not currently implemented");
3837   ierr = (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr);
3838   PetscFunctionReturn(0);
3839 }
3840 
3841 #undef __FUNC__
3842 #define __FUNC__ /*<a name=""></a>*/"MatGetMaps"
3843 /*@C
3844    MatGetMaps - Returns the maps associated with the matrix.
3845 
3846    Not Collective
3847 
3848    Input Parameter:
3849 .  mat - the matrix
3850 
3851    Output Parameters:
3852 +  rmap - the row (right) map
3853 -  cmap - the column (left) map
3854 
3855    Level: developer
3856 
3857 .keywords: matrix, get, map
3858 @*/
3859 int MatGetMaps(Mat mat,Map *rmap,Map *cmap)
3860 {
3861   int ierr;
3862 
3863   PetscFunctionBegin;
3864   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3865   ierr = (*mat->ops->getmaps)(mat,rmap,cmap);CHKERRQ(ierr);
3866   PetscFunctionReturn(0);
3867 }
3868 
3869 /*
3870       Version that works for all PETSc matrices
3871 */
3872 #undef __FUNC__
3873 #define __FUNC__ /*<a name=""></a>*/"MatGetMaps_Petsc"
3874 int MatGetMaps_Petsc(Mat mat,Map *rmap,Map *cmap)
3875 {
3876   PetscFunctionBegin;
3877   if (rmap) *rmap = mat->rmap;
3878   if (cmap) *cmap = mat->cmap;
3879   PetscFunctionReturn(0);
3880 }
3881 
3882 #undef __FUNC__
3883 #define __FUNC__ /*<a name=""></a>*/"MatSetStashInitialSize"
3884 /*@
3885    MatSetStashInitialSize - sets the sizes of the matrix stash, that is
3886    used during the assembly process to store values that belong to
3887    other processors.
3888 
3889    Not Collective
3890 
3891    Input Parameters:
3892 +  mat   - the matrix
3893 .  size  - the initial size of the stash.
3894 -  bsize - the initial size of the block-stash(if used).
3895 
3896    Options Database Keys:
3897 +   -matstash_initial_size <size> or <size0,size1,...sizep-1>
3898 -   -matstash_block_initial_size <bsize>  or <bsize0,bsize1,...bsizep-1>
3899 
3900    Level: intermediate
3901 
3902    Notes:
3903      The block-stash is used for values set with VecSetValuesBlocked() while
3904      the stash is used for values set with VecSetValues()
3905 
3906      Run with the option -log_info and look for output of the form
3907      MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
3908      to determine the appropriate value, MM, to use for size and
3909      MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
3910      to determine the value, BMM to use for bsize
3911 
3912 .keywords: matrix, stash, assembly
3913 @*/
3914 int MatSetStashInitialSize(Mat mat,int size, int bsize)
3915 {
3916   int ierr;
3917 
3918   PetscFunctionBegin;
3919   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3920   ierr = MatStashSetInitialSize_Private(&mat->stash,size);CHKERRQ(ierr);
3921   ierr = MatStashSetInitialSize_Private(&mat->bstash,bsize);CHKERRQ(ierr);
3922   PetscFunctionReturn(0);
3923 }
3924 
3925 #undef __FUNC__
3926 #define __FUNC__ /*<a name=""></a>*/"MatInterpolateAdd"
3927 /*@
3928    MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of
3929      the matrix
3930 
3931    Collective on Mat
3932 
3933    Input Parameters:
3934 +  mat   - the matrix
3935 .  x,y - the vectors
3936 -  w - where the result is stored
3937 
3938    Level: intermediate
3939 
3940    Notes:
3941     w may be the same vector as y.
3942 
3943     This allows one to use either the restriction or interpolation (its transpose)
3944     matrix to do the interpolation
3945 
3946 .keywords: interpolate,
3947 
3948 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
3949 
3950 @*/
3951 int MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w)
3952 {
3953   int M,N,ierr;
3954 
3955   PetscFunctionBegin;
3956   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
3957   if (N > M) {
3958     ierr = MatMultTransposeAdd(A,x,y,w);CHKERRQ(ierr);
3959   } else {
3960     ierr = MatMultAdd(A,x,y,w);CHKERRQ(ierr);
3961   }
3962   PetscFunctionReturn(0);
3963 }
3964 
3965 #undef __FUNC__
3966 #define __FUNC__ /*<a name=""></a>*/"MatInterpolate"
3967 /*@
3968    MatInterpolate - y = A*x or A'*x depending on the shape of
3969      the matrix
3970 
3971    Collective on Mat
3972 
3973    Input Parameters:
3974 +  mat   - the matrix
3975 -  x,y - the vectors
3976 
3977    Level: intermediate
3978 
3979    Notes:
3980     This allows one to use either the restriction or interpolation (its transpose)
3981     matrix to do the interpolation
3982 
3983 .keywords: interpolate,
3984 
3985 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
3986 
3987 @*/
3988 int MatInterpolate(Mat A,Vec x,Vec y)
3989 {
3990   int M,N,ierr;
3991 
3992   PetscFunctionBegin;
3993   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
3994   if (N > M) {
3995     ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
3996   } else {
3997     ierr = MatMult(A,x,y);CHKERRQ(ierr);
3998   }
3999   PetscFunctionReturn(0);
4000 }
4001 
4002 #undef __FUNC__
4003 #define __FUNC__ /*<a name=""></a>*/"MatRestrict"
4004 /*@
4005    MatRestrict - y = A*x or A'*x
4006 
4007    Collective on Mat
4008 
4009    Input Parameters:
4010 +  mat   - the matrix
4011 -  x,y - the vectors
4012 
4013    Level: intermediate
4014 
4015    Notes:
4016     This allows one to use either the restriction or interpolation (its transpose)
4017     matrix to do the restriction
4018 
4019 .keywords: interpolate,
4020 
4021 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatInterpolate()
4022 
4023 @*/
4024 int MatRestrict(Mat A,Vec x,Vec y)
4025 {
4026   int M,N,ierr;
4027 
4028   PetscFunctionBegin;
4029   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
4030   if (N > M) {
4031     ierr = MatMult(A,x,y);CHKERRQ(ierr);
4032   } else {
4033     ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
4034   }
4035   PetscFunctionReturn(0);
4036 }
4037