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