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