xref: /petsc/src/mat/interface/matrix.c (revision da1d93b3f862e692ef185b13ee9cc501acd81181)
1 /*$Id: matrix.c,v 1.368 2000/05/16 17:12:43 bsmith Exp bsmith $*/
2 
3 /*
4    This is where the abstract matrix operations are defined
5 */
6 
7 #include "src/mat/matimpl.h"        /*I "petscmat.h" I*/
8 #include "src/vec/vecimpl.h"
9 
10 #undef __FUNC__
11 #define __FUNC__ /*<a name=""></a>*/"MatGetRow"
12 /*@C
13    MatGetRow - Gets a row of a matrix.  You MUST call MatRestoreRow()
14    for each row that you get to ensure that your application does
15    not bleed memory.
16 
17    Not Collective
18 
19    Input Parameters:
20 +  mat - the matrix
21 -  row - the row to get
22 
23    Output Parameters:
24 +  ncols -  the number of nonzeros in the row
25 .  cols - if 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    Level: intermediate
2753 
2754 .keywords: matrix, option, row-oriented, column-oriented, sorted, nonzero
2755 @*/
2756 int MatSetOption(Mat mat,MatOption op)
2757 {
2758   int ierr;
2759 
2760   PetscFunctionBegin;
2761   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2762   if (mat->ops->setoption) {ierr = (*mat->ops->setoption)(mat,op);CHKERRQ(ierr);}
2763   PetscFunctionReturn(0);
2764 }
2765 
2766 #undef __FUNC__
2767 #define __FUNC__ /*<a name=""></a>*/"MatZeroEntries"
2768 /*@
2769    MatZeroEntries - Zeros all entries of a matrix.  For sparse matrices
2770    this routine retains the old nonzero structure.
2771 
2772    Collective on Mat
2773 
2774    Input Parameters:
2775 .  mat - the matrix
2776 
2777    Level: intermediate
2778 
2779 .keywords: matrix, zero, entries
2780 
2781 .seealso: MatZeroRows()
2782 @*/
2783 int MatZeroEntries(Mat mat)
2784 {
2785   int ierr;
2786 
2787   PetscFunctionBegin;
2788   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2789   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2790   if (!mat->ops->zeroentries) SETERRQ(PETSC_ERR_SUP,0,"");
2791 
2792   PLogEventBegin(MAT_ZeroEntries,mat,0,0,0);
2793   ierr = (*mat->ops->zeroentries)(mat);CHKERRQ(ierr);
2794   PLogEventEnd(MAT_ZeroEntries,mat,0,0,0);
2795   PetscFunctionReturn(0);
2796 }
2797 
2798 #undef __FUNC__
2799 #define __FUNC__ /*<a name=""></a>*/"MatZeroRows"
2800 /*@C
2801    MatZeroRows - Zeros all entries (except possibly the main diagonal)
2802    of a set of rows of a matrix.
2803 
2804    Collective on Mat
2805 
2806    Input Parameters:
2807 +  mat - the matrix
2808 .  is - index set of rows to remove
2809 -  diag - pointer to value put in all diagonals of eliminated rows.
2810           Note that diag is not a pointer to an array, but merely a
2811           pointer to a single value.
2812 
2813    Notes:
2814    For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
2815    but does not release memory.  For the dense and block diagonal
2816    formats this does not alter the nonzero structure.
2817 
2818    If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
2819    of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
2820    merely zeroed.
2821 
2822    The user can set a value in the diagonal entry (or for the AIJ and
2823    row formats can optionally remove the main diagonal entry from the
2824    nonzero structure as well, by passing a null pointer (PETSC_NULL
2825    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
2826 
2827    For the parallel case, all processes that share the matrix (i.e.,
2828    those in the communicator used for matrix creation) MUST call this
2829    routine, regardless of whether any rows being zeroed are owned by
2830    them.
2831 
2832 
2833    Level: intermediate
2834 
2835 .keywords: matrix, zero, rows, boundary conditions
2836 
2837 .seealso: MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
2838 @*/
2839 int MatZeroRows(Mat mat,IS is,Scalar *diag)
2840 {
2841   int ierr;
2842 
2843   PetscFunctionBegin;
2844   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2845   PetscValidHeaderSpecific(is,IS_COOKIE);
2846   if (diag) PetscValidScalarPointer(diag);
2847   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2848   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2849   if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,0,"");
2850 
2851   ierr = (*mat->ops->zerorows)(mat,is,diag);CHKERRQ(ierr);
2852   ierr = MatView_Private(mat);CHKERRQ(ierr);
2853   PetscFunctionReturn(0);
2854 }
2855 
2856 #undef __FUNC__
2857 #define __FUNC__ /*<a name=""></a>*/"MatZeroRowsLocal"
2858 /*@C
2859    MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
2860    of a set of rows of a matrix; using local numbering of rows.
2861 
2862    Collective on Mat
2863 
2864    Input Parameters:
2865 +  mat - the matrix
2866 .  is - index set of rows to remove
2867 -  diag - pointer to value put in all diagonals of eliminated rows.
2868           Note that diag is not a pointer to an array, but merely a
2869           pointer to a single value.
2870 
2871    Notes:
2872    Before calling MatZeroRowsLocal(), the user must first set the
2873    local-to-global mapping by calling MatSetLocalToGlobalMapping().
2874 
2875    For the AIJ matrix formats this removes the old nonzero structure,
2876    but does not release memory.  For the dense and block diagonal
2877    formats this does not alter the nonzero structure.
2878 
2879    The user can set a value in the diagonal entry (or for the AIJ and
2880    row formats can optionally remove the main diagonal entry from the
2881    nonzero structure as well, by passing a null pointer (PETSC_NULL
2882    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
2883 
2884    Level: intermediate
2885 
2886 .keywords: matrix, zero, rows, boundary conditions
2887 
2888 .seealso: MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
2889 @*/
2890 int MatZeroRowsLocal(Mat mat,IS is,Scalar *diag)
2891 {
2892   int ierr;
2893   IS  newis;
2894 
2895   PetscFunctionBegin;
2896   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2897   PetscValidHeaderSpecific(is,IS_COOKIE);
2898   if (diag) PetscValidScalarPointer(diag);
2899   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2900   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2901   if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,0,"");
2902   if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Need to provide local to global mapping to matrix first");
2903 
2904   ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);CHKERRQ(ierr);
2905   ierr = (*mat->ops->zerorows)(mat,newis,diag);CHKERRQ(ierr);
2906   ierr = ISDestroy(newis);CHKERRQ(ierr);
2907   PetscFunctionReturn(0);
2908 }
2909 
2910 #undef __FUNC__
2911 #define __FUNC__ /*<a name=""></a>*/"MatGetSize"
2912 /*@
2913    MatGetSize - Returns the numbers of rows and columns in a matrix.
2914 
2915    Not Collective
2916 
2917    Input Parameter:
2918 .  mat - the matrix
2919 
2920    Output Parameters:
2921 +  m - the number of global rows
2922 -  n - the number of global columns
2923 
2924    Level: beginner
2925 
2926 .keywords: matrix, dimension, size, rows, columns, global, get
2927 
2928 .seealso: MatGetLocalSize()
2929 @*/
2930 int MatGetSize(Mat mat,int *m,int* n)
2931 {
2932   int ierr;
2933 
2934   PetscFunctionBegin;
2935   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2936   ierr = (*mat->ops->getsize)(mat,m,n);CHKERRQ(ierr);
2937   PetscFunctionReturn(0);
2938 }
2939 
2940 #undef __FUNC__
2941 #define __FUNC__ /*<a name=""></a>*/"MatGetLocalSize"
2942 /*@
2943    MatGetLocalSize - Returns the number of rows and columns in a matrix
2944    stored locally.  This information may be implementation dependent, so
2945    use with care.
2946 
2947    Not Collective
2948 
2949    Input Parameters:
2950 .  mat - the matrix
2951 
2952    Output Parameters:
2953 +  m - the number of local rows
2954 -  n - the number of local columns
2955 
2956    Level: beginner
2957 
2958 .keywords: matrix, dimension, size, local, rows, columns, get
2959 
2960 .seealso: MatGetSize()
2961 @*/
2962 int MatGetLocalSize(Mat mat,int *m,int* n)
2963 {
2964   int ierr;
2965 
2966   PetscFunctionBegin;
2967   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2968   ierr = (*mat->ops->getlocalsize)(mat,m,n);CHKERRQ(ierr);
2969   PetscFunctionReturn(0);
2970 }
2971 
2972 #undef __FUNC__
2973 #define __FUNC__ /*<a name=""></a>*/"MatGetOwnershipRange"
2974 /*@
2975    MatGetOwnershipRange - Returns the range of matrix rows owned by
2976    this processor, assuming that the matrix is laid out with the first
2977    n1 rows on the first processor, the next n2 rows on the second, etc.
2978    For certain parallel layouts this range may not be well defined.
2979 
2980    Not Collective
2981 
2982    Input Parameters:
2983 .  mat - the matrix
2984 
2985    Output Parameters:
2986 +  m - the global index of the first local row
2987 -  n - one more than the global index of the last local row
2988 
2989    Level: beginner
2990 
2991 .keywords: matrix, get, range, ownership
2992 @*/
2993 int MatGetOwnershipRange(Mat mat,int *m,int* n)
2994 {
2995   int ierr;
2996 
2997   PetscFunctionBegin;
2998   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2999   if (m) PetscValidIntPointer(m);
3000   if (n) PetscValidIntPointer(n);
3001   if (!mat->ops->getownershiprange) SETERRQ(PETSC_ERR_SUP,0,"");
3002   ierr = (*mat->ops->getownershiprange)(mat,m,n);CHKERRQ(ierr);
3003   PetscFunctionReturn(0);
3004 }
3005 
3006 #undef __FUNC__
3007 #define __FUNC__ /*<a name=""></a>*/"MatILUFactorSymbolic"
3008 /*@
3009    MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
3010    Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric()
3011    to complete the factorization.
3012 
3013    Collective on Mat
3014 
3015    Input Parameters:
3016 +  mat - the matrix
3017 .  row - row permutation
3018 .  column - column permutation
3019 -  info - structure containing
3020 $      levels - number of levels of fill.
3021 $      expected fill - as ratio of original fill.
3022 $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
3023                 missing diagonal entries)
3024 
3025    Output Parameters:
3026 .  fact - new matrix that has been symbolically factored
3027 
3028    Notes:
3029    See the users manual for additional information about
3030    choosing the fill factor for better efficiency.
3031 
3032    Most users should employ the simplified SLES interface for linear solvers
3033    instead of working directly with matrix algebra routines such as this.
3034    See, e.g., SLESCreate().
3035 
3036    Level: developer
3037 
3038 .keywords: matrix, factor, incomplete, ILU, symbolic, fill
3039 
3040 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
3041           MatGetOrdering()
3042 
3043 @*/
3044 int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatILUInfo *info,Mat *fact)
3045 {
3046   int ierr;
3047 
3048   PetscFunctionBegin;
3049   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3050   PetscValidPointer(fact);
3051   if (info && info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Levels of fill negative %d",info->levels);
3052   if (info && info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",info->fill);
3053   if (!mat->ops->ilufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Only MatCreateMPIRowbs() matrices support parallel ILU");
3054   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3055   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3056 
3057   PLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);
3058   ierr = (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr);
3059   PLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);
3060   PetscFunctionReturn(0);
3061 }
3062 
3063 #undef __FUNC__
3064 #define __FUNC__ /*<a name=""></a>*/"MatIncompleteCholeskyFactorSymbolic"
3065 /*@
3066    MatIncompleteCholeskyFactorSymbolic - Performs symbolic incomplete
3067    Cholesky factorization for a symmetric matrix.  Use
3068    MatCholeskyFactorNumeric() to complete the factorization.
3069 
3070    Collective on Mat
3071 
3072    Input Parameters:
3073 +  mat - the matrix
3074 .  perm - row and column permutation
3075 .  fill - levels of fill
3076 -  f - expected fill as ratio of original fill
3077 
3078    Output Parameter:
3079 .  fact - the factored matrix
3080 
3081    Notes:
3082    Currently only no-fill factorization is supported.
3083 
3084    Most users should employ the simplified SLES interface for linear solvers
3085    instead of working directly with matrix algebra routines such as this.
3086    See, e.g., SLESCreate().
3087 
3088    Level: developer
3089 
3090 .keywords: matrix, factor, incomplete, ICC, Cholesky, symbolic, fill
3091 
3092 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor()
3093 @*/
3094 int MatIncompleteCholeskyFactorSymbolic(Mat mat,IS perm,PetscReal f,int fill,Mat *fact)
3095 {
3096   int ierr;
3097 
3098   PetscFunctionBegin;
3099   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3100   PetscValidPointer(fact);
3101   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3102   if (fill < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Fill negative %d",fill);
3103   if (f < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",f);
3104   if (!mat->ops->incompletecholeskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support ICC in parallel");
3105   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3106 
3107   PLogEventBegin(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);
3108   ierr = (*mat->ops->incompletecholeskyfactorsymbolic)(mat,perm,f,fill,fact);CHKERRQ(ierr);
3109   PLogEventEnd(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);
3110   PetscFunctionReturn(0);
3111 }
3112 
3113 #undef __FUNC__
3114 #define __FUNC__ /*<a name=""></a>*/"MatGetArray"
3115 /*@C
3116    MatGetArray - Returns a pointer to the element values in the matrix.
3117    The result of this routine is dependent on the underlying matrix data
3118    structure, and may not even work for certain matrix types.  You MUST
3119    call MatRestoreArray() when you no longer need to access the array.
3120 
3121    Not Collective
3122 
3123    Input Parameter:
3124 .  mat - the matrix
3125 
3126    Output Parameter:
3127 .  v - the location of the values
3128 
3129    Currently returns an array only for the dense formats, giving access to
3130    the local portion of the matrix in the usual Fortran column-oriented format.
3131 
3132    Fortran Note:
3133    This routine is used differently from Fortran, e.g.,
3134 .vb
3135         Mat         mat
3136         Scalar      mat_array(1)
3137         PetscOffset i_mat
3138         int         ierr
3139         call MatGetArray(mat,mat_array,i_mat,ierr)
3140 
3141   C  Access first local entry in matrix; note that array is
3142   C  treated as one dimensional
3143         value = mat_array(i_mat + 1)
3144 
3145         [... other code ...]
3146         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3147 .ve
3148 
3149    See the Fortran chapter of the users manual and
3150    petsc/src/mat/examples/tests for details.
3151 
3152    Level: advanced
3153 
3154 .keywords: matrix, array, elements, values
3155 
3156 .seealso: MatRestoreArray(), MatGetArrayF90()
3157 @*/
3158 int MatGetArray(Mat mat,Scalar **v)
3159 {
3160   int ierr;
3161 
3162   PetscFunctionBegin;
3163   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3164   PetscValidPointer(v);
3165   if (!mat->ops->getarray) SETERRQ(PETSC_ERR_SUP,0,"");
3166   ierr = (*mat->ops->getarray)(mat,v);CHKERRQ(ierr);
3167   PetscFunctionReturn(0);
3168 }
3169 
3170 #undef __FUNC__
3171 #define __FUNC__ /*<a name=""></a>*/"MatRestoreArray"
3172 /*@C
3173    MatRestoreArray - Restores the matrix after MatGetArray() has been called.
3174 
3175    Not Collective
3176 
3177    Input Parameter:
3178 +  mat - the matrix
3179 -  v - the location of the values
3180 
3181    Fortran Note:
3182    This routine is used differently from Fortran, e.g.,
3183 .vb
3184         Mat         mat
3185         Scalar      mat_array(1)
3186         PetscOffset i_mat
3187         int         ierr
3188         call MatGetArray(mat,mat_array,i_mat,ierr)
3189 
3190   C  Access first local entry in matrix; note that array is
3191   C  treated as one dimensional
3192         value = mat_array(i_mat + 1)
3193 
3194         [... other code ...]
3195         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3196 .ve
3197 
3198    See the Fortran chapter of the users manual and
3199    petsc/src/mat/examples/tests for details
3200 
3201    Level: advanced
3202 
3203 .keywords: matrix, array, elements, values, restore
3204 
3205 .seealso: MatGetArray(), MatRestoreArrayF90()
3206 @*/
3207 int MatRestoreArray(Mat mat,Scalar **v)
3208 {
3209   int ierr;
3210 
3211   PetscFunctionBegin;
3212   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3213   PetscValidPointer(v);
3214   if (!mat->ops->restorearray) SETERRQ(PETSC_ERR_SUP,0,"");
3215   ierr = (*mat->ops->restorearray)(mat,v);CHKERRQ(ierr);
3216   PetscFunctionReturn(0);
3217 }
3218 
3219 #undef __FUNC__
3220 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrices"
3221 /*@C
3222    MatGetSubMatrices - Extracts several submatrices from a matrix. If submat
3223    points to an array of valid matrices, they may be reused to store the new
3224    submatrices.
3225 
3226    Collective on Mat
3227 
3228    Input Parameters:
3229 +  mat - the matrix
3230 .  n   - the number of submatrixes to be extracted (on this processor, may be zero)
3231 .  irow, icol - index sets of rows and columns to extract
3232 -  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3233 
3234    Output Parameter:
3235 .  submat - the array of submatrices
3236 
3237    Notes:
3238    MatGetSubMatrices() can extract only sequential submatrices
3239    (from both sequential and parallel matrices). Use MatGetSubMatrix()
3240    to extract a parallel submatrix.
3241 
3242    When extracting submatrices from a parallel matrix, each processor can
3243    form a different submatrix by setting the rows and columns of its
3244    individual index sets according to the local submatrix desired.
3245 
3246    When finished using the submatrices, the user should destroy
3247    them with MatDestroySubMatrices().
3248 
3249    MAT_REUSE_MATRIX can only be used when the nonzero structure of the
3250    original matrix has not changed from that last call to MatGetSubMatrices().
3251 
3252    Fortran Note:
3253    The Fortran interface is slightly different from that given below; it
3254    requires one to pass in  as submat a Mat (integer) array of size at least m.
3255 
3256    Level: advanced
3257 
3258 .keywords: matrix, get, submatrix, submatrices
3259 
3260 .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal()
3261 @*/
3262 int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatReuse scall,Mat **submat)
3263 {
3264   int ierr;
3265 
3266   PetscFunctionBegin;
3267   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3268   if (!mat->ops->getsubmatrices) SETERRQ(PETSC_ERR_SUP,0,"");
3269   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3270 
3271   PLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);
3272   ierr = (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);CHKERRQ(ierr);
3273   PLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);
3274 
3275   PetscFunctionReturn(0);
3276 }
3277 
3278 #undef __FUNC__
3279 #define __FUNC__ /*<a name=""></a>*/"MatDestroyMatrices"
3280 /*@C
3281    MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices().
3282 
3283    Collective on Mat
3284 
3285    Input Parameters:
3286 +  n - the number of local matrices
3287 -  mat - the matrices
3288 
3289    Level: advanced
3290 
3291 .keywords: matrix, destroy, submatrix, submatrices
3292 
3293 .seealso: MatGetSubMatrices()
3294 @*/
3295 int MatDestroyMatrices(int n,Mat **mat)
3296 {
3297   int ierr,i;
3298 
3299   PetscFunctionBegin;
3300   if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,1,"Trying to destroy negative number of matrices %d",n);
3301   PetscValidPointer(mat);
3302   for (i=0; i<n; i++) {
3303     ierr = MatDestroy((*mat)[i]);CHKERRQ(ierr);
3304   }
3305   if (n) {ierr = PetscFree(*mat);CHKERRQ(ierr);}
3306   PetscFunctionReturn(0);
3307 }
3308 
3309 #undef __FUNC__
3310 #define __FUNC__ /*<a name=""></a>*/"MatIncreaseOverlap"
3311 /*@
3312    MatIncreaseOverlap - Given a set of submatrices indicated by index sets,
3313    replaces the index sets by larger ones that represent submatrices with
3314    additional overlap.
3315 
3316    Collective on Mat
3317 
3318    Input Parameters:
3319 +  mat - the matrix
3320 .  n   - the number of index sets
3321 .  is  - the array of pointers to index sets
3322 -  ov  - the additional overlap requested
3323 
3324    Level: developer
3325 
3326 .keywords: matrix, overlap, Schwarz
3327 
3328 .seealso: MatGetSubMatrices()
3329 @*/
3330 int MatIncreaseOverlap(Mat mat,int n,IS *is,int ov)
3331 {
3332   int ierr;
3333 
3334   PetscFunctionBegin;
3335   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3336   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3337   if (mat->factor)     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3338 
3339   if (!ov) PetscFunctionReturn(0);
3340   if (!mat->ops->increaseoverlap) SETERRQ(PETSC_ERR_SUP,0,"");
3341   PLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);
3342   ierr = (*mat->ops->increaseoverlap)(mat,n,is,ov);CHKERRQ(ierr);
3343   PLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);
3344   PetscFunctionReturn(0);
3345 }
3346 
3347 #undef __FUNC__
3348 #define __FUNC__ /*<a name=""></a>*/"MatPrintHelp"
3349 /*@
3350    MatPrintHelp - Prints all the options for the matrix.
3351 
3352    Collective on Mat
3353 
3354    Input Parameter:
3355 .  mat - the matrix
3356 
3357    Options Database Keys:
3358 +  -help - Prints matrix options
3359 -  -h - Prints matrix options
3360 
3361    Level: developer
3362 
3363 .keywords: mat, help
3364 
3365 .seealso: MatCreate(), MatCreateXXX()
3366 @*/
3367 int MatPrintHelp(Mat mat)
3368 {
3369   static PetscTruth called = PETSC_FALSE;
3370   int               ierr;
3371   MPI_Comm          comm;
3372 
3373   PetscFunctionBegin;
3374   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3375 
3376   comm = mat->comm;
3377   if (!called) {
3378     ierr = (*PetscHelpPrintf)(comm,"General matrix options:\n");CHKERRQ(ierr);
3379     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_info: view basic matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr);
3380     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_info_detailed: view detailed matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr);
3381     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_draw: draw nonzero matrix structure during MatAssemblyEnd()\n");CHKERRQ(ierr);
3382     ierr = (*PetscHelpPrintf)(comm,"      -draw_pause <sec>: set seconds of display pause\n");CHKERRQ(ierr);
3383     ierr = (*PetscHelpPrintf)(comm,"      -display <name>: set alternate display\n");CHKERRQ(ierr);
3384     called = PETSC_TRUE;
3385   }
3386   if (mat->ops->printhelp) {
3387     ierr = (*mat->ops->printhelp)(mat);CHKERRQ(ierr);
3388   }
3389   PetscFunctionReturn(0);
3390 }
3391 
3392 #undef __FUNC__
3393 #define __FUNC__ /*<a name=""></a>*/"MatGetBlockSize"
3394 /*@
3395    MatGetBlockSize - Returns the matrix block size; useful especially for the
3396    block row and block diagonal formats.
3397 
3398    Not Collective
3399 
3400    Input Parameter:
3401 .  mat - the matrix
3402 
3403    Output Parameter:
3404 .  bs - block size
3405 
3406    Notes:
3407    Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG.
3408    Block row formats are MATSEQBAIJ, MATMPIBAIJ
3409 
3410    Level: intermediate
3411 
3412 .keywords: matrix, get, block, size
3413 
3414 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag()
3415 @*/
3416 int MatGetBlockSize(Mat mat,int *bs)
3417 {
3418   int ierr;
3419 
3420   PetscFunctionBegin;
3421   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3422   PetscValidIntPointer(bs);
3423   if (!mat->ops->getblocksize) SETERRQ(PETSC_ERR_SUP,0,"");
3424   ierr = (*mat->ops->getblocksize)(mat,bs);CHKERRQ(ierr);
3425   PetscFunctionReturn(0);
3426 }
3427 
3428 #undef __FUNC__
3429 #define __FUNC__ /*<a name=""></a>*/"MatGetRowIJ"
3430 /*@C
3431     MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices.
3432 
3433    Collective on Mat
3434 
3435     Input Parameters:
3436 +   mat - the matrix
3437 .   shift -  0 or 1 indicating we want the indices starting at 0 or 1
3438 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3439                 symmetrized
3440 
3441     Output Parameters:
3442 +   n - number of rows in the (possibly compressed) matrix
3443 .   ia - the row pointers
3444 .   ja - the column indices
3445 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3446 
3447     Level: developer
3448 
3449 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3450 @*/
3451 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3452 {
3453   int ierr;
3454 
3455   PetscFunctionBegin;
3456   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3457   if (ia) PetscValidIntPointer(ia);
3458   if (ja) PetscValidIntPointer(ja);
3459   PetscValidIntPointer(done);
3460   if (!mat->ops->getrowij) *done = PETSC_FALSE;
3461   else {
3462     *done = PETSC_TRUE;
3463     ierr  = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3464   }
3465   PetscFunctionReturn(0);
3466 }
3467 
3468 #undef __FUNC__
3469 #define __FUNC__ /*<a name=""></a>*/"MatGetColumnIJ"
3470 /*@C
3471     MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices.
3472 
3473     Collective on Mat
3474 
3475     Input Parameters:
3476 +   mat - the matrix
3477 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3478 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3479                 symmetrized
3480 
3481     Output Parameters:
3482 +   n - number of columns in the (possibly compressed) matrix
3483 .   ia - the column pointers
3484 .   ja - the row indices
3485 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3486 
3487     Level: developer
3488 
3489 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3490 @*/
3491 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3492 {
3493   int ierr;
3494 
3495   PetscFunctionBegin;
3496   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3497   if (ia) PetscValidIntPointer(ia);
3498   if (ja) PetscValidIntPointer(ja);
3499   PetscValidIntPointer(done);
3500 
3501   if (!mat->ops->getcolumnij) *done = PETSC_FALSE;
3502   else {
3503     *done = PETSC_TRUE;
3504     ierr  = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3505   }
3506   PetscFunctionReturn(0);
3507 }
3508 
3509 #undef __FUNC__
3510 #define __FUNC__ /*<a name=""></a>*/"MatRestoreRowIJ"
3511 /*@C
3512     MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with
3513     MatGetRowIJ().
3514 
3515     Collective on Mat
3516 
3517     Input Parameters:
3518 +   mat - the matrix
3519 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3520 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3521                 symmetrized
3522 
3523     Output Parameters:
3524 +   n - size of (possibly compressed) matrix
3525 .   ia - the row pointers
3526 .   ja - the column indices
3527 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3528 
3529     Level: developer
3530 
3531 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3532 @*/
3533 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3534 {
3535   int ierr;
3536 
3537   PetscFunctionBegin;
3538   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3539   if (ia) PetscValidIntPointer(ia);
3540   if (ja) PetscValidIntPointer(ja);
3541   PetscValidIntPointer(done);
3542 
3543   if (!mat->ops->restorerowij) *done = PETSC_FALSE;
3544   else {
3545     *done = PETSC_TRUE;
3546     ierr  = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3547   }
3548   PetscFunctionReturn(0);
3549 }
3550 
3551 #undef __FUNC__
3552 #define __FUNC__ /*<a name=""></a>*/"MatRestoreColumnIJ"
3553 /*@C
3554     MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with
3555     MatGetColumnIJ().
3556 
3557     Collective on Mat
3558 
3559     Input Parameters:
3560 +   mat - the matrix
3561 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3562 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3563                 symmetrized
3564 
3565     Output Parameters:
3566 +   n - size of (possibly compressed) matrix
3567 .   ia - the column pointers
3568 .   ja - the row indices
3569 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3570 
3571     Level: developer
3572 
3573 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3574 @*/
3575 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3576 {
3577   int ierr;
3578 
3579   PetscFunctionBegin;
3580   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3581   if (ia) PetscValidIntPointer(ia);
3582   if (ja) PetscValidIntPointer(ja);
3583   PetscValidIntPointer(done);
3584 
3585   if (!mat->ops->restorecolumnij) *done = PETSC_FALSE;
3586   else {
3587     *done = PETSC_TRUE;
3588     ierr  = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3589   }
3590   PetscFunctionReturn(0);
3591 }
3592 
3593 #undef __FUNC__
3594 #define __FUNC__ /*<a name=""></a>*/"MatColoringPatch"
3595 /*@C
3596     MatColoringPatch -Used inside matrix coloring routines that
3597     use MatGetRowIJ() and/or MatGetColumnIJ().
3598 
3599     Collective on Mat
3600 
3601     Input Parameters:
3602 +   mat - the matrix
3603 .   n   - number of colors
3604 -   colorarray - array indicating color for each column
3605 
3606     Output Parameters:
3607 .   iscoloring - coloring generated using colorarray information
3608 
3609     Level: developer
3610 
3611 .seealso: MatGetRowIJ(), MatGetColumnIJ()
3612 
3613 .keywords: mat, coloring, patch
3614 @*/
3615 int MatColoringPatch(Mat mat,int n,int *colorarray,ISColoring *iscoloring)
3616 {
3617   int ierr;
3618 
3619   PetscFunctionBegin;
3620   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3621   PetscValidIntPointer(colorarray);
3622 
3623   if (!mat->ops->coloringpatch) {SETERRQ(PETSC_ERR_SUP,0,"");}
3624   else {
3625     ierr  = (*mat->ops->coloringpatch)(mat,n,colorarray,iscoloring);CHKERRQ(ierr);
3626   }
3627   PetscFunctionReturn(0);
3628 }
3629 
3630 
3631 #undef __FUNC__
3632 #define __FUNC__ /*<a name=""></a>*/"MatSetUnfactored"
3633 /*@
3634    MatSetUnfactored - Resets a factored matrix to be treated as unfactored.
3635 
3636    Collective on Mat
3637 
3638    Input Parameter:
3639 .  mat - the factored matrix to be reset
3640 
3641    Notes:
3642    This routine should be used only with factored matrices formed by in-place
3643    factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE
3644    format).  This option can save memory, for example, when solving nonlinear
3645    systems with a matrix-free Newton-Krylov method and a matrix-based, in-place
3646    ILU(0) preconditioner.
3647 
3648    Note that one can specify in-place ILU(0) factorization by calling
3649 .vb
3650      PCType(pc,PCILU);
3651      PCILUSeUseInPlace(pc);
3652 .ve
3653    or by using the options -pc_type ilu -pc_ilu_in_place
3654 
3655    In-place factorization ILU(0) can also be used as a local
3656    solver for the blocks within the block Jacobi or additive Schwarz
3657    methods (runtime option: -sub_pc_ilu_in_place).  See the discussion
3658    of these preconditioners in the users manual for details on setting
3659    local solver options.
3660 
3661    Most users should employ the simplified SLES interface for linear solvers
3662    instead of working directly with matrix algebra routines such as this.
3663    See, e.g., SLESCreate().
3664 
3665    Level: developer
3666 
3667 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace()
3668 
3669 .keywords: matrix-free, in-place ILU, in-place LU
3670 @*/
3671 int MatSetUnfactored(Mat mat)
3672 {
3673   int ierr;
3674 
3675   PetscFunctionBegin;
3676   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3677   mat->factor = 0;
3678   if (!mat->ops->setunfactored) PetscFunctionReturn(0);
3679   ierr = (*mat->ops->setunfactored)(mat);CHKERRQ(ierr);
3680   PetscFunctionReturn(0);
3681 }
3682 
3683 #undef __FUNC__
3684 #define __FUNC__ /*<a name=""></a>*/"MatGetType"
3685 /*@C
3686    MatGetType - Gets the matrix type and name (as a string) from the matrix.
3687 
3688    Not Collective
3689 
3690    Input Parameter:
3691 .  mat - the matrix
3692 
3693    Output Parameter:
3694 +  type - the matrix type (or use PETSC_NULL)
3695 -  name - name of matrix type (or use PETSC_NULL)
3696 
3697    Level: intermediate
3698 
3699 .keywords: matrix, get, type, name
3700 @*/
3701 int MatGetType(Mat mat,MatType *type,char **name)
3702 {
3703   int  itype = (int)mat->type;
3704   char *matname[10];
3705 
3706   PetscFunctionBegin;
3707   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3708 
3709   if (type) *type = (MatType) mat->type;
3710   if (name) {
3711     /* Note:  Be sure that this list corresponds to the enum in petscmat.h */
3712     matname[0] = "MATSEQDENSE";
3713     matname[1] = "MATSEQAIJ";
3714     matname[2] = "MATMPIAIJ";
3715     matname[3] = "MATSHELL";
3716     matname[4] = "MATMPIROWBS";
3717     matname[5] = "MATSEQBDIAG";
3718     matname[6] = "MATMPIBDIAG";
3719     matname[7] = "MATMPIDENSE";
3720     matname[8] = "MATSEQBAIJ";
3721     matname[9] = "MATMPIBAIJ";
3722 
3723     if (itype < 0 || itype > 9) *name = "Unknown matrix type";
3724     else                        *name = matname[itype];
3725   }
3726   PetscFunctionReturn(0);
3727 }
3728 
3729 /*MC
3730     MatGetArrayF90 - Accesses a matrix array from Fortran90.
3731 
3732     Synopsis:
3733     MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3734 
3735     Not collective
3736 
3737     Input Parameter:
3738 .   x - matrix
3739 
3740     Output Parameters:
3741 +   xx_v - the Fortran90 pointer to the array
3742 -   ierr - error code
3743 
3744     Example of Usage:
3745 .vb
3746       Scalar, pointer xx_v(:)
3747       ....
3748       call MatGetArrayF90(x,xx_v,ierr)
3749       a = xx_v(3)
3750       call MatRestoreArrayF90(x,xx_v,ierr)
3751 .ve
3752 
3753     Notes:
3754     Not yet supported for all F90 compilers
3755 
3756     Level: advanced
3757 
3758 .seealso:  MatRestoreArrayF90(), MatGetArray(), MatRestoreArray()
3759 
3760 .keywords:  matrix, array, f90
3761 M*/
3762 
3763 /*MC
3764     MatRestoreArrayF90 - Restores a matrix array that has been
3765     accessed with MatGetArrayF90().
3766 
3767     Synopsis:
3768     MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3769 
3770     Not collective
3771 
3772     Input Parameters:
3773 +   x - matrix
3774 -   xx_v - the Fortran90 pointer to the array
3775 
3776     Output Parameter:
3777 .   ierr - error code
3778 
3779     Example of Usage:
3780 .vb
3781        Scalar, pointer xx_v(:)
3782        ....
3783        call MatGetArrayF90(x,xx_v,ierr)
3784        a = xx_v(3)
3785        call MatRestoreArrayF90(x,xx_v,ierr)
3786 .ve
3787 
3788     Notes:
3789     Not yet supported for all F90 compilers
3790 
3791     Level: advanced
3792 
3793 .seealso:  MatGetArrayF90(), MatGetArray(), MatRestoreArray()
3794 
3795 .keywords:  matrix, array, f90
3796 M*/
3797 
3798 
3799 #undef __FUNC__
3800 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrix"
3801 /*@
3802     MatGetSubMatrix - Gets a single submatrix on the same number of processors
3803                       as the original matrix.
3804 
3805     Collective on Mat
3806 
3807     Input Parameters:
3808 +   mat - the original matrix
3809 .   isrow - rows this processor should obtain
3810 .   iscol - columns for all processors you wish to keep
3811 .   csize - number of columns "local" to this processor (does nothing for sequential
3812             matrices). This should match the result from VecGetLocalSize(x,...) if you
3813             plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE
3814 -   cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3815 
3816     Output Parameter:
3817 .   newmat - the new submatrix, of the same type as the old
3818 
3819     Level: advanced
3820 
3821 .keywords: matrix, get, submatrix, submatrices
3822 
3823 .seealso: MatGetSubMatrices()
3824 @*/
3825 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat)
3826 {
3827   int     ierr, size;
3828   Mat     *local;
3829 
3830   PetscFunctionBegin;
3831   ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr);
3832 
3833   /* if original matrix is on just one processor then use submatrix generated */
3834   if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) {
3835     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr);
3836     PetscFunctionReturn(0);
3837   } else if (!mat->ops->getsubmatrix && size == 1) {
3838     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr);
3839     *newmat = *local;
3840     ierr = PetscFree(local);CHKERRQ(ierr);
3841     PetscFunctionReturn(0);
3842   }
3843 
3844   if (!mat->ops->getsubmatrix) SETERRQ(PETSC_ERR_SUP,0,"Not currently implemented");
3845   ierr = (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr);
3846   PetscFunctionReturn(0);
3847 }
3848 
3849 #undef __FUNC__
3850 #define __FUNC__ /*<a name=""></a>*/"MatGetMaps"
3851 /*@C
3852    MatGetMaps - Returns the maps associated with the matrix.
3853 
3854    Not Collective
3855 
3856    Input Parameter:
3857 .  mat - the matrix
3858 
3859    Output Parameters:
3860 +  rmap - the row (right) map
3861 -  cmap - the column (left) map
3862 
3863    Level: developer
3864 
3865 .keywords: matrix, get, map
3866 @*/
3867 int MatGetMaps(Mat mat,Map *rmap,Map *cmap)
3868 {
3869   int ierr;
3870 
3871   PetscFunctionBegin;
3872   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3873   ierr = (*mat->ops->getmaps)(mat,rmap,cmap);CHKERRQ(ierr);
3874   PetscFunctionReturn(0);
3875 }
3876 
3877 /*
3878       Version that works for all PETSc matrices
3879 */
3880 #undef __FUNC__
3881 #define __FUNC__ /*<a name=""></a>*/"MatGetMaps_Petsc"
3882 int MatGetMaps_Petsc(Mat mat,Map *rmap,Map *cmap)
3883 {
3884   PetscFunctionBegin;
3885   if (rmap) *rmap = mat->rmap;
3886   if (cmap) *cmap = mat->cmap;
3887   PetscFunctionReturn(0);
3888 }
3889 
3890 #undef __FUNC__
3891 #define __FUNC__ /*<a name=""></a>*/"MatSetStashInitialSize"
3892 /*@
3893    MatSetStashInitialSize - sets the sizes of the matrix stash, that is
3894    used during the assembly process to store values that belong to
3895    other processors.
3896 
3897    Not Collective
3898 
3899    Input Parameters:
3900 +  mat   - the matrix
3901 .  size  - the initial size of the stash.
3902 -  bsize - the initial size of the block-stash(if used).
3903 
3904    Options Database Keys:
3905 +   -matstash_initial_size <size> or <size0,size1,...sizep-1>
3906 -   -matstash_block_initial_size <bsize>  or <bsize0,bsize1,...bsizep-1>
3907 
3908    Level: intermediate
3909 
3910    Notes:
3911      The block-stash is used for values set with VecSetValuesBlocked() while
3912      the stash is used for values set with VecSetValues()
3913 
3914      Run with the option -log_info and look for output of the form
3915      MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
3916      to determine the appropriate value, MM, to use for size and
3917      MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
3918      to determine the value, BMM to use for bsize
3919 
3920 .keywords: matrix, stash, assembly
3921 @*/
3922 int MatSetStashInitialSize(Mat mat,int size, int bsize)
3923 {
3924   int ierr;
3925 
3926   PetscFunctionBegin;
3927   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3928   ierr = MatStashSetInitialSize_Private(&mat->stash,size);CHKERRQ(ierr);
3929   ierr = MatStashSetInitialSize_Private(&mat->bstash,bsize);CHKERRQ(ierr);
3930   PetscFunctionReturn(0);
3931 }
3932 
3933 #undef __FUNC__
3934 #define __FUNC__ /*<a name=""></a>*/"MatInterpolateAdd"
3935 /*@
3936    MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of
3937      the matrix
3938 
3939    Collective on Mat
3940 
3941    Input Parameters:
3942 +  mat   - the matrix
3943 .  x,y - the vectors
3944 -  w - where the result is stored
3945 
3946    Level: intermediate
3947 
3948    Notes:
3949     w may be the same vector as y.
3950 
3951     This allows one to use either the restriction or interpolation (its transpose)
3952     matrix to do the interpolation
3953 
3954 .keywords: interpolate,
3955 
3956 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
3957 
3958 @*/
3959 int MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w)
3960 {
3961   int M,N,ierr;
3962 
3963   PetscFunctionBegin;
3964   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
3965   if (N > M) {
3966     ierr = MatMultTransposeAdd(A,x,y,w);CHKERRQ(ierr);
3967   } else {
3968     ierr = MatMultAdd(A,x,y,w);CHKERRQ(ierr);
3969   }
3970   PetscFunctionReturn(0);
3971 }
3972 
3973 #undef __FUNC__
3974 #define __FUNC__ /*<a name=""></a>*/"MatInterpolate"
3975 /*@
3976    MatInterpolate - y = A*x or A'*x depending on the shape of
3977      the matrix
3978 
3979    Collective on Mat
3980 
3981    Input Parameters:
3982 +  mat   - the matrix
3983 -  x,y - the vectors
3984 
3985    Level: intermediate
3986 
3987    Notes:
3988     This allows one to use either the restriction or interpolation (its transpose)
3989     matrix to do the interpolation
3990 
3991 .keywords: interpolate,
3992 
3993 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
3994 
3995 @*/
3996 int MatInterpolate(Mat A,Vec x,Vec y)
3997 {
3998   int M,N,ierr;
3999 
4000   PetscFunctionBegin;
4001   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
4002   if (N > M) {
4003     ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
4004   } else {
4005     ierr = MatMult(A,x,y);CHKERRQ(ierr);
4006   }
4007   PetscFunctionReturn(0);
4008 }
4009 
4010 #undef __FUNC__
4011 #define __FUNC__ /*<a name=""></a>*/"MatRestrict"
4012 /*@
4013    MatRestrict - y = A*x or A'*x
4014 
4015    Collective on Mat
4016 
4017    Input Parameters:
4018 +  mat   - the matrix
4019 -  x,y - the vectors
4020 
4021    Level: intermediate
4022 
4023    Notes:
4024     This allows one to use either the restriction or interpolation (its transpose)
4025     matrix to do the restriction
4026 
4027 .keywords: interpolate,
4028 
4029 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatInterpolate()
4030 
4031 @*/
4032 int MatRestrict(Mat A,Vec x,Vec y)
4033 {
4034   int M,N,ierr;
4035 
4036   PetscFunctionBegin;
4037   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
4038   if (N > M) {
4039     ierr = MatMult(A,x,y);CHKERRQ(ierr);
4040   } else {
4041     ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
4042   }
4043   PetscFunctionReturn(0);
4044 }
4045