xref: /petsc/src/mat/interface/matrix.c (revision 5d1a10b1e93e00413872aeb7c1118d31f04cd00a)
1 /*$Id: matrix.c,v 1.375 2000/06/06 19:01:58 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   if (x->ops->setlocaltoglobalmapping) {
679     ierr = (*x->ops->setlocaltoglobalmapping)(x,mapping);CHKERRQ(ierr);
680   } else {
681     x->mapping = mapping;
682     ierr = PetscObjectReference((PetscObject)mapping);CHKERRQ(ierr);
683   }
684   PetscFunctionReturn(0);
685 }
686 
687 #undef __FUNC__
688 #define __FUNC__ /*<a name=""></a>*/"MatSetLocalToGlobalMappingBlock"
689 /*@
690    MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use
691    by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
692    entries using a local (per-processor) numbering.
693 
694    Not Collective
695 
696    Input Parameters:
697 +  x - the matrix
698 -  mapping - mapping created with ISLocalToGlobalMappingCreate() or
699              ISLocalToGlobalMappingCreateIS()
700 
701    Level: intermediate
702 
703 .keywords: matrix, set, values, local ordering
704 
705 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
706            MatSetValuesBlocked(), MatSetValuesLocal()
707 @*/
708 int MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping mapping)
709 {
710   int ierr;
711   PetscFunctionBegin;
712   PetscValidHeaderSpecific(x,MAT_COOKIE);
713   PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE);
714   if (x->bmapping) {
715     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Mapping already set for matrix");
716   }
717 
718   x->bmapping = mapping;
719   ierr = PetscObjectReference((PetscObject)mapping);CHKERRQ(ierr);
720   PetscFunctionReturn(0);
721 }
722 
723 #undef __FUNC__
724 #define __FUNC__ /*<a name=""></a>*/"MatSetValuesLocal"
725 /*@
726    MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
727    using a local ordering of the nodes.
728 
729    Not Collective
730 
731    Input Parameters:
732 +  x - the matrix
733 .  nrow, irow - number of rows and their local indices
734 .  ncol, icol - number of columns and their local indices
735 .  y -  a logically two-dimensional array of values
736 -  addv - either INSERT_VALUES or ADD_VALUES, where
737    ADD_VALUES adds values to any existing entries, and
738    INSERT_VALUES replaces existing entries with new values
739 
740    Notes:
741    Before calling MatSetValuesLocal(), the user must first set the
742    local-to-global mapping by calling MatSetLocalToGlobalMapping().
743 
744    Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES
745    options cannot be mixed without intervening calls to the assembly
746    routines.
747 
748    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
749    MUST be called after all calls to MatSetValuesLocal() have been completed.
750 
751    Level: intermediate
752 
753 .keywords: matrix, set, values, local ordering
754 
755 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping()
756 @*/
757 int MatSetValuesLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,Scalar *y,InsertMode addv)
758 {
759   int ierr,irowm[2048],icolm[2048];
760 
761   PetscFunctionBegin;
762   PetscValidHeaderSpecific(mat,MAT_COOKIE);
763   PetscValidIntPointer(irow);
764   PetscValidIntPointer(icol);
765   PetscValidScalarPointer(y);
766 
767   if (mat->insertmode == NOT_SET_VALUES) {
768     mat->insertmode = addv;
769   }
770 #if defined(PETSC_USE_BOPT_g)
771   else if (mat->insertmode != addv) {
772     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values");
773   }
774   if (!mat->ops->setvalueslocal && (nrow > 2048 || ncol > 2048)) {
775     SETERRQ2(PETSC_ERR_SUP,0,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
776   }
777   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
778 #endif
779 
780   if (mat->assembled) {
781     mat->was_assembled = PETSC_TRUE;
782     mat->assembled     = PETSC_FALSE;
783   }
784   PLogEventBegin(MAT_SetValues,mat,0,0,0);
785   if (!mat->ops->setvalueslocal) {
786     ierr = ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm);CHKERRQ(ierr);
787     ierr = ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);CHKERRQ(ierr);
788     ierr = (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr);
789   } else {
790     ierr = (*mat->ops->setvalueslocal)(mat,nrow,irow,ncol,icol,y,addv);CHKERRQ(ierr);
791   }
792   PLogEventEnd(MAT_SetValues,mat,0,0,0);
793   PetscFunctionReturn(0);
794 }
795 
796 #undef __FUNC__
797 #define __FUNC__ /*<a name=""></a>*/"MatSetValuesBlockedLocal"
798 /*@
799    MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
800    using a local ordering of the nodes a block at a time.
801 
802    Not Collective
803 
804    Input Parameters:
805 +  x - the matrix
806 .  nrow, irow - number of rows and their local indices
807 .  ncol, icol - number of columns and their local indices
808 .  y -  a logically two-dimensional array of values
809 -  addv - either INSERT_VALUES or ADD_VALUES, where
810    ADD_VALUES adds values to any existing entries, and
811    INSERT_VALUES replaces existing entries with new values
812 
813    Notes:
814    Before calling MatSetValuesBlockedLocal(), the user must first set the
815    local-to-global mapping by calling MatSetLocalToGlobalMappingBlock(),
816    where the mapping MUST be set for matrix blocks, not for matrix elements.
817 
818    Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
819    options cannot be mixed without intervening calls to the assembly
820    routines.
821 
822    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
823    MUST be called after all calls to MatSetValuesBlockedLocal() have been completed.
824 
825    Level: intermediate
826 
827 .keywords: matrix, set, values, blocked, local
828 
829 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked()
830 @*/
831 int MatSetValuesBlockedLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,Scalar *y,InsertMode addv)
832 {
833   int ierr,irowm[2048],icolm[2048];
834 
835   PetscFunctionBegin;
836   PetscValidHeaderSpecific(mat,MAT_COOKIE);
837   PetscValidIntPointer(irow);
838   PetscValidIntPointer(icol);
839   PetscValidScalarPointer(y);
840   if (mat->insertmode == NOT_SET_VALUES) {
841     mat->insertmode = addv;
842   }
843 #if defined(PETSC_USE_BOPT_g)
844   else if (mat->insertmode != addv) {
845     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values");
846   }
847   if (!mat->bmapping) {
848     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Local to global never set with MatSetLocalToGlobalMappingBlock()");
849   }
850   if (nrow > 2048 || ncol > 2048) {
851     SETERRQ2(PETSC_ERR_SUP,0,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
852   }
853   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
854 #endif
855 
856   if (mat->assembled) {
857     mat->was_assembled = PETSC_TRUE;
858     mat->assembled     = PETSC_FALSE;
859   }
860   PLogEventBegin(MAT_SetValues,mat,0,0,0);
861   ierr = ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm);CHKERRQ(ierr);
862   ierr = ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm);CHKERRQ(ierr);
863   ierr = (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr);
864   PLogEventEnd(MAT_SetValues,mat,0,0,0);
865   PetscFunctionReturn(0);
866 }
867 
868 /* --------------------------------------------------------*/
869 #undef __FUNC__
870 #define __FUNC__ /*<a name=""></a>*/"MatMult"
871 /*@
872    MatMult - Computes the matrix-vector product, y = Ax.
873 
874    Collective on Mat and Vec
875 
876    Input Parameters:
877 +  mat - the matrix
878 -  x   - the vector to be multilplied
879 
880    Output Parameters:
881 .  y - the result
882 
883    Notes:
884    The vectors x and y cannot be the same.  I.e., one cannot
885    call MatMult(A,y,y).
886 
887    Level: beginner
888 
889 .keywords: matrix, multiply, matrix-vector product
890 
891 .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
892 @*/
893 int MatMult(Mat mat,Vec x,Vec y)
894 {
895   int ierr;
896 
897   PetscFunctionBegin;
898   PetscValidHeaderSpecific(mat,MAT_COOKIE);
899   PetscValidHeaderSpecific(x,VEC_COOKIE);
900   PetscValidHeaderSpecific(y,VEC_COOKIE);
901 
902   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
903   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
904   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"x and y must be different vectors");
905   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
906   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
907   if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: local dim %d %d",mat->m,y->n);
908 
909   PLogEventBegin(MAT_Mult,mat,x,y,0);
910   ierr = (*mat->ops->mult)(mat,x,y);CHKERRQ(ierr);
911   PLogEventEnd(MAT_Mult,mat,x,y,0);
912 
913   PetscFunctionReturn(0);
914 }
915 
916 #undef __FUNC__
917 #define __FUNC__ /*<a name=""></a>*/"MatMultTranspose"
918 /*@
919    MatMultTranspose - Computes matrix transpose times a vector.
920 
921    Collective on Mat and Vec
922 
923    Input Parameters:
924 +  mat - the matrix
925 -  x   - the vector to be multilplied
926 
927    Output Parameters:
928 .  y - the result
929 
930    Notes:
931    The vectors x and y cannot be the same.  I.e., one cannot
932    call MatMultTranspose(A,y,y).
933 
934    Level: beginner
935 
936 .keywords: matrix, multiply, matrix-vector product, transpose
937 
938 .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd()
939 @*/
940 int MatMultTranspose(Mat mat,Vec x,Vec y)
941 {
942   int ierr;
943 
944   PetscFunctionBegin;
945   PetscValidHeaderSpecific(mat,MAT_COOKIE);
946   PetscValidHeaderSpecific(x,VEC_COOKIE);
947   PetscValidHeaderSpecific(y,VEC_COOKIE);
948 
949   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
950   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
951   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"x and y must be different vectors");
952   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
953   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
954 
955   PLogEventBegin(MAT_MultTranspose,mat,x,y,0);
956   ierr = (*mat->ops->multtranspose)(mat,x,y);CHKERRQ(ierr);
957   PLogEventEnd(MAT_MultTranspose,mat,x,y,0);
958   PetscFunctionReturn(0);
959 }
960 
961 #undef __FUNC__
962 #define __FUNC__ /*<a name=""></a>*/"MatMultAdd"
963 /*@
964     MatMultAdd -  Computes v3 = v2 + A * v1.
965 
966     Collective on Mat and Vec
967 
968     Input Parameters:
969 +   mat - the matrix
970 -   v1, v2 - the vectors
971 
972     Output Parameters:
973 .   v3 - the result
974 
975     Notes:
976     The vectors v1 and v3 cannot be the same.  I.e., one cannot
977     call MatMultAdd(A,v1,v2,v1).
978 
979     Level: beginner
980 
981 .keywords: matrix, multiply, matrix-vector product, add
982 
983 .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd()
984 @*/
985 int MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
986 {
987   int ierr;
988 
989   PetscFunctionBegin;
990   PetscValidHeaderSpecific(mat,MAT_COOKIE);
991   PetscValidHeaderSpecific(v1,VEC_COOKIE);
992   PetscValidHeaderSpecific(v2,VEC_COOKIE);
993   PetscValidHeaderSpecific(v3,VEC_COOKIE);
994 
995   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
996   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
997   if (mat->N != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v1: global dim %d %d",mat->N,v1->N);
998   if (mat->M != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: global dim %d %d",mat->M,v2->N);
999   if (mat->M != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: global dim %d %d",mat->M,v3->N);
1000   if (mat->m != v3->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: local dim %d %d",mat->m,v3->n);
1001   if (mat->m != v2->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: local dim %d %d",mat->m,v2->n);
1002   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,0,"v1 and v3 must be different vectors");
1003 
1004   PLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);
1005   ierr = (*mat->ops->multadd)(mat,v1,v2,v3);CHKERRQ(ierr);
1006   PLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);
1007   PetscFunctionReturn(0);
1008 }
1009 
1010 #undef __FUNC__
1011 #define __FUNC__ /*<a name=""></a>*/"MatMultTransposeAdd"
1012 /*@
1013    MatMultTransposeAdd - Computes v3 = v2 + A' * v1.
1014 
1015    Collective on Mat and Vec
1016 
1017    Input Parameters:
1018 +  mat - the matrix
1019 -  v1, v2 - the vectors
1020 
1021    Output Parameters:
1022 .  v3 - the result
1023 
1024    Notes:
1025    The vectors v1 and v3 cannot be the same.  I.e., one cannot
1026    call MatMultTransposeAdd(A,v1,v2,v1).
1027 
1028    Level: beginner
1029 
1030 .keywords: matrix, multiply, matrix-vector product, transpose, add
1031 
1032 .seealso: MatMultTranspose(), MatMultAdd(), MatMult()
1033 @*/
1034 int MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1035 {
1036   int ierr;
1037 
1038   PetscFunctionBegin;
1039   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1040   PetscValidHeaderSpecific(v1,VEC_COOKIE);
1041   PetscValidHeaderSpecific(v2,VEC_COOKIE);
1042   PetscValidHeaderSpecific(v3,VEC_COOKIE);
1043 
1044   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1045   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1046   if (!mat->ops->multtransposeadd) SETERRQ(PETSC_ERR_SUP,0,"");
1047   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,0,"v1 and v3 must be different vectors");
1048   if (mat->M != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v1: global dim %d %d",mat->M,v1->N);
1049   if (mat->N != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: global dim %d %d",mat->N,v2->N);
1050   if (mat->N != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: global dim %d %d",mat->N,v3->N);
1051 
1052   PLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);
1053   ierr = (*mat->ops->multtransposeadd)(mat,v1,v2,v3);CHKERRQ(ierr);
1054   PLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);
1055   PetscFunctionReturn(0);
1056 }
1057 /* ------------------------------------------------------------*/
1058 #undef __FUNC__
1059 #define __FUNC__ /*<a name=""></a>*/"MatGetInfo"
1060 /*@C
1061    MatGetInfo - Returns information about matrix storage (number of
1062    nonzeros, memory, etc.).
1063 
1064    Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used
1065    as the flag
1066 
1067    Input Parameters:
1068 .  mat - the matrix
1069 
1070    Output Parameters:
1071 +  flag - flag indicating the type of parameters to be returned
1072    (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
1073    MAT_GLOBAL_SUM - sum over all processors)
1074 -  info - matrix information context
1075 
1076    Notes:
1077    The MatInfo context contains a variety of matrix data, including
1078    number of nonzeros allocated and used, number of mallocs during
1079    matrix assembly, etc.  Additional information for factored matrices
1080    is provided (such as the fill ratio, number of mallocs during
1081    factorization, etc.).  Much of this info is printed to STDOUT
1082    when using the runtime options
1083 $       -log_info -mat_view_info
1084 
1085    Example for C/C++ Users:
1086    See the file ${PETSC_DIR}/include/petscmat.h for a complete list of
1087    data within the MatInfo context.  For example,
1088 .vb
1089       MatInfo info;
1090       Mat     A;
1091       double  mal, nz_a, nz_u;
1092 
1093       MatGetInfo(A,MAT_LOCAL,&info);
1094       mal  = info.mallocs;
1095       nz_a = info.nz_allocated;
1096 .ve
1097 
1098    Example for Fortran Users:
1099    Fortran users should declare info as a double precision
1100    array of dimension MAT_INFO_SIZE, and then extract the parameters
1101    of interest.  See the file ${PETSC_DIR}/include/finclude/petscmat.h
1102    a complete list of parameter names.
1103 .vb
1104       double  precision info(MAT_INFO_SIZE)
1105       double  precision mal, nz_a
1106       Mat     A
1107       integer ierr
1108 
1109       call MatGetInfo(A,MAT_LOCAL,info,ierr)
1110       mal = info(MAT_INFO_MALLOCS)
1111       nz_a = info(MAT_INFO_NZ_ALLOCATED)
1112 .ve
1113 
1114     Level: intermediate
1115 
1116 .keywords: matrix, get, info, storage, nonzeros, memory, fill
1117 @*/
1118 int MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
1119 {
1120   int ierr;
1121 
1122   PetscFunctionBegin;
1123   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1124   PetscValidPointer(info);
1125   if (!mat->ops->getinfo) SETERRQ(PETSC_ERR_SUP,0,"");
1126   ierr = (*mat->ops->getinfo)(mat,flag,info);CHKERRQ(ierr);
1127   PetscFunctionReturn(0);
1128 }
1129 
1130 /* ----------------------------------------------------------*/
1131 #undef __FUNC__
1132 #define __FUNC__ /*<a name=""></a>*/"MatILUDTFactor"
1133 /*@C
1134    MatILUDTFactor - Performs a drop tolerance ILU factorization.
1135 
1136    Collective on Mat
1137 
1138    Input Parameters:
1139 +  mat - the matrix
1140 .  info - information about the factorization to be done
1141 .  row - row permutation
1142 -  col - column permutation
1143 
1144    Output Parameters:
1145 .  fact - the factored matrix
1146 
1147    Level: developer
1148 
1149    Notes:
1150    Most users should employ the simplified SLES interface for linear solvers
1151    instead of working directly with matrix algebra routines such as this.
1152    See, e.g., SLESCreate().
1153 
1154    This is currently only supported for the SeqAIJ matrix format using code
1155    from Yousef Saad's SPARSEKIT2  package (translated to C with f2c) and/or
1156    Matlab. SPARSEKIT2 is copyrighted by Yousef Saad with the GNU copyright
1157    and thus can be distributed with PETSc.
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   ierr = OptionsHasName(mat->prefix,"-mat_view_binary",&flg);CHKERRQ(ierr);
2560   if (flg) {
2561     ierr = MatView(mat,VIEWER_BINARY_(mat->comm));CHKERRQ(ierr);
2562     ierr = ViewerFlush(VIEWER_BINARY_(mat->comm));CHKERRQ(ierr);
2563   }
2564   PetscFunctionReturn(0);
2565 }
2566 
2567 #undef __FUNC__
2568 #define __FUNC__ /*<a name=""></a>*/"MatAssemblyEnd"
2569 /*@
2570    MatAssemblyEnd - Completes assembling the matrix.  This routine should
2571    be called after MatAssemblyBegin().
2572 
2573    Collective on Mat
2574 
2575    Input Parameters:
2576 +  mat - the matrix
2577 -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
2578 
2579    Options Database Keys:
2580 +  -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
2581 .  -mat_view_info_detailed - Prints more detailed info
2582 .  -mat_view - Prints matrix in ASCII format
2583 .  -mat_view_matlab - Prints matrix in Matlab format
2584 .  -mat_view_draw - Draws nonzero structure of matrix, using MatView() and DrawOpenX().
2585 .  -display <name> - Sets display name (default is host)
2586 -  -draw_pause <sec> - Sets number of seconds to pause after display
2587 
2588    Notes:
2589    MatSetValues() generally caches the values.  The matrix is ready to
2590    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
2591    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
2592    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
2593    using the matrix.
2594 
2595    Level: beginner
2596 
2597 .keywords: matrix, assembly, assemble, end
2598 
2599 .seealso: MatAssemblyBegin(), MatSetValues(), DrawOpenX(), MatView(), MatAssembled()
2600 @*/
2601 int MatAssemblyEnd(Mat mat,MatAssemblyType type)
2602 {
2603   int        ierr;
2604   static int inassm = 0;
2605 
2606   PetscFunctionBegin;
2607   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2608 
2609   inassm++;
2610   MatAssemblyEnd_InUse++;
2611   if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
2612     PLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);
2613     if (mat->ops->assemblyend) {
2614       ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr);
2615     }
2616     PLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);
2617   } else {
2618     if (mat->ops->assemblyend) {
2619       ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr);
2620     }
2621   }
2622 
2623   /* Flush assembly is not a true assembly */
2624   if (type != MAT_FLUSH_ASSEMBLY) {
2625     mat->assembled  = PETSC_TRUE; mat->num_ass++;
2626   }
2627   mat->insertmode = NOT_SET_VALUES;
2628   MatAssemblyEnd_InUse--;
2629 
2630   if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) {
2631     ierr = MatView_Private(mat);CHKERRQ(ierr);
2632   }
2633   inassm--;
2634   PetscFunctionReturn(0);
2635 }
2636 
2637 
2638 #undef __FUNC__
2639 #define __FUNC__ /*<a name=""></a>*/"MatCompress"
2640 /*@
2641    MatCompress - Tries to store the matrix in as little space as
2642    possible.  May fail if memory is already fully used, since it
2643    tries to allocate new space.
2644 
2645    Collective on Mat
2646 
2647    Input Parameters:
2648 .  mat - the matrix
2649 
2650    Level: advanced
2651 
2652 .keywords: matrix, compress
2653 @*/
2654 int MatCompress(Mat mat)
2655 {
2656   int ierr;
2657 
2658   PetscFunctionBegin;
2659   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2660   if (mat->ops->compress) {ierr = (*mat->ops->compress)(mat);CHKERRQ(ierr);}
2661   PetscFunctionReturn(0);
2662 }
2663 
2664 #undef __FUNC__
2665 #define __FUNC__ /*<a name=""></a>*/"MatSetOption"
2666 /*@
2667    MatSetOption - Sets a parameter option for a matrix. Some options
2668    may be specific to certain storage formats.  Some options
2669    determine how values will be inserted (or added). Sorted,
2670    row-oriented input will generally assemble the fastest. The default
2671    is row-oriented, nonsorted input.
2672 
2673    Collective on Mat
2674 
2675    Input Parameters:
2676 +  mat - the matrix
2677 -  option - the option, one of those listed below (and possibly others),
2678              e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR
2679 
2680    Options Describing Matrix Structure:
2681 +    MAT_SYMMETRIC - symmetric in terms of both structure and value
2682 -    MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure
2683 
2684    Options For Use with MatSetValues():
2685    Insert a logically dense subblock, which can be
2686 +    MAT_ROW_ORIENTED - row-oriented
2687 .    MAT_COLUMN_ORIENTED - column-oriented
2688 .    MAT_ROWS_SORTED - sorted by row
2689 .    MAT_ROWS_UNSORTED - not sorted by row
2690 .    MAT_COLUMNS_SORTED - sorted by column
2691 -    MAT_COLUMNS_UNSORTED - not sorted by column
2692 
2693    Not these options reflect the data you pass in with MatSetValues(); it has
2694    nothing to do with how the data is stored internally in the matrix
2695    data structure.
2696 
2697    When (re)assembling a matrix, we can restrict the input for
2698    efficiency/debugging purposes.  These options include
2699 +    MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be
2700         allowed if they generate a new nonzero
2701 .    MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed
2702 .    MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if
2703          they generate a nonzero in a new diagonal (for block diagonal format only)
2704 .    MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only)
2705 .    MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries
2706 .    MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry
2707 -    MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly
2708 
2709    Notes:
2710    Some options are relevant only for particular matrix types and
2711    are thus ignored by others.  Other options are not supported by
2712    certain matrix types and will generate an error message if set.
2713 
2714    If using a Fortran 77 module to compute a matrix, one may need to
2715    use the column-oriented option (or convert to the row-oriented
2716    format).
2717 
2718    MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion
2719    that would generate a new entry in the nonzero structure is instead
2720    ignored.  Thus, if memory has not alredy been allocated for this particular
2721    data, then the insertion is ignored. For dense matrices, in which
2722    the entire array is allocated, no entries are ever ignored.
2723 
2724    MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion
2725    that would generate a new entry in the nonzero structure instead produces
2726    an error. (Currently supported for AIJ and BAIJ formats only.)
2727    This is a useful flag when using SAME_NONZERO_PATTERN in calling
2728    SLESSetOperators() to ensure that the nonzero pattern truely does
2729    remain unchanged.
2730 
2731    MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion
2732    that would generate a new entry that has not been preallocated will
2733    instead produce an error. (Currently supported for AIJ and BAIJ formats
2734    only.) This is a useful flag when debugging matrix memory preallocation.
2735 
2736    MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for
2737    other processors should be dropped, rather than stashed.
2738    This is useful if you know that the "owning" processor is also
2739    always generating the correct matrix entries, so that PETSc need
2740    not transfer duplicate entries generated on another processor.
2741 
2742    MAT_USE_HASH_TABLE indicates that a hash table be used to improve the
2743    searches during matrix assembly. When this flag is set, the hash table
2744    is created during the first Matrix Assembly. This hash table is
2745    used the next time through, during MatSetVaules()/MatSetVaulesBlocked()
2746    to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag
2747    should be used with MAT_USE_HASH_TABLE flag. This option is currently
2748    supported by MATMPIBAIJ format only.
2749 
2750    MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries
2751    are kept in the nonzero structure
2752 
2753    MAT_IGNORE_ZERO_ENTRIES - when using ADD_VALUES for AIJ matrices this will stop
2754    zero values from creating a zero location in the matrix
2755 
2756    MAT_USE_INODES - indicates using inode version of the code - works with AIJ and
2757    ROWBS matrix types
2758 
2759    MAT_DO_NOT_USE_INODES - indicates not using inode version of the code - works
2760    with AIJ and ROWBS matrix types
2761 
2762    Level: intermediate
2763 
2764 .keywords: matrix, option, row-oriented, column-oriented, sorted, nonzero
2765 @*/
2766 int MatSetOption(Mat mat,MatOption op)
2767 {
2768   int ierr;
2769 
2770   PetscFunctionBegin;
2771   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2772   if (mat->ops->setoption) {ierr = (*mat->ops->setoption)(mat,op);CHKERRQ(ierr);}
2773   PetscFunctionReturn(0);
2774 }
2775 
2776 #undef __FUNC__
2777 #define __FUNC__ /*<a name=""></a>*/"MatZeroEntries"
2778 /*@
2779    MatZeroEntries - Zeros all entries of a matrix.  For sparse matrices
2780    this routine retains the old nonzero structure.
2781 
2782    Collective on Mat
2783 
2784    Input Parameters:
2785 .  mat - the matrix
2786 
2787    Level: intermediate
2788 
2789 .keywords: matrix, zero, entries
2790 
2791 .seealso: MatZeroRows()
2792 @*/
2793 int MatZeroEntries(Mat mat)
2794 {
2795   int ierr;
2796 
2797   PetscFunctionBegin;
2798   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2799   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2800   if (!mat->ops->zeroentries) SETERRQ(PETSC_ERR_SUP,0,"");
2801 
2802   PLogEventBegin(MAT_ZeroEntries,mat,0,0,0);
2803   ierr = (*mat->ops->zeroentries)(mat);CHKERRQ(ierr);
2804   PLogEventEnd(MAT_ZeroEntries,mat,0,0,0);
2805   PetscFunctionReturn(0);
2806 }
2807 
2808 #undef __FUNC__
2809 #define __FUNC__ /*<a name=""></a>*/"MatZeroRows"
2810 /*@C
2811    MatZeroRows - Zeros all entries (except possibly the main diagonal)
2812    of a set of rows of a matrix.
2813 
2814    Collective on Mat
2815 
2816    Input Parameters:
2817 +  mat - the matrix
2818 .  is - index set of rows to remove
2819 -  diag - pointer to value put in all diagonals of eliminated rows.
2820           Note that diag is not a pointer to an array, but merely a
2821           pointer to a single value.
2822 
2823    Notes:
2824    For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
2825    but does not release memory.  For the dense and block diagonal
2826    formats this does not alter the nonzero structure.
2827 
2828    If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
2829    of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
2830    merely zeroed.
2831 
2832    The user can set a value in the diagonal entry (or for the AIJ and
2833    row formats can optionally remove the main diagonal entry from the
2834    nonzero structure as well, by passing a null pointer (PETSC_NULL
2835    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
2836 
2837    For the parallel case, all processes that share the matrix (i.e.,
2838    those in the communicator used for matrix creation) MUST call this
2839    routine, regardless of whether any rows being zeroed are owned by
2840    them.
2841 
2842 
2843    Level: intermediate
2844 
2845 .keywords: matrix, zero, rows, boundary conditions
2846 
2847 .seealso: MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
2848 @*/
2849 int MatZeroRows(Mat mat,IS is,Scalar *diag)
2850 {
2851   int ierr;
2852 
2853   PetscFunctionBegin;
2854   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2855   PetscValidHeaderSpecific(is,IS_COOKIE);
2856   if (diag) PetscValidScalarPointer(diag);
2857   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2858   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2859   if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,0,"");
2860 
2861   ierr = (*mat->ops->zerorows)(mat,is,diag);CHKERRQ(ierr);
2862   ierr = MatView_Private(mat);CHKERRQ(ierr);
2863   PetscFunctionReturn(0);
2864 }
2865 
2866 #undef __FUNC__
2867 #define __FUNC__ /*<a name=""></a>*/"MatZeroRowsLocal"
2868 /*@C
2869    MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
2870    of a set of rows of a matrix; using local numbering of rows.
2871 
2872    Collective on Mat
2873 
2874    Input Parameters:
2875 +  mat - the matrix
2876 .  is - index set of rows to remove
2877 -  diag - pointer to value put in all diagonals of eliminated rows.
2878           Note that diag is not a pointer to an array, but merely a
2879           pointer to a single value.
2880 
2881    Notes:
2882    Before calling MatZeroRowsLocal(), the user must first set the
2883    local-to-global mapping by calling MatSetLocalToGlobalMapping().
2884 
2885    For the AIJ matrix formats this removes the old nonzero structure,
2886    but does not release memory.  For the dense and block diagonal
2887    formats this does not alter the nonzero structure.
2888 
2889    The user can set a value in the diagonal entry (or for the AIJ and
2890    row formats can optionally remove the main diagonal entry from the
2891    nonzero structure as well, by passing a null pointer (PETSC_NULL
2892    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
2893 
2894    Level: intermediate
2895 
2896 .keywords: matrix, zero, rows, boundary conditions
2897 
2898 .seealso: MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
2899 @*/
2900 int MatZeroRowsLocal(Mat mat,IS is,Scalar *diag)
2901 {
2902   int ierr;
2903   IS  newis;
2904 
2905   PetscFunctionBegin;
2906   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2907   PetscValidHeaderSpecific(is,IS_COOKIE);
2908   if (diag) PetscValidScalarPointer(diag);
2909   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2910   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2911 
2912   if (mat->ops->zerorowslocal) {
2913     ierr = (*mat->ops->zerorowslocal)(mat,is,diag);CHKERRQ(ierr);
2914   } else {
2915     if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Need to provide local to global mapping to matrix first");
2916     ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);CHKERRQ(ierr);
2917     ierr = (*mat->ops->zerorows)(mat,newis,diag);CHKERRQ(ierr);
2918     ierr = ISDestroy(newis);CHKERRQ(ierr);
2919   }
2920   PetscFunctionReturn(0);
2921 }
2922 
2923 #undef __FUNC__
2924 #define __FUNC__ /*<a name=""></a>*/"MatGetSize"
2925 /*@
2926    MatGetSize - Returns the numbers of rows and columns in a matrix.
2927 
2928    Not Collective
2929 
2930    Input Parameter:
2931 .  mat - the matrix
2932 
2933    Output Parameters:
2934 +  m - the number of global rows
2935 -  n - the number of global columns
2936 
2937    Level: beginner
2938 
2939 .keywords: matrix, dimension, size, rows, columns, global, get
2940 
2941 .seealso: MatGetLocalSize()
2942 @*/
2943 int MatGetSize(Mat mat,int *m,int* n)
2944 {
2945   int ierr;
2946 
2947   PetscFunctionBegin;
2948   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2949   ierr = (*mat->ops->getsize)(mat,m,n);CHKERRQ(ierr);
2950   PetscFunctionReturn(0);
2951 }
2952 
2953 #undef __FUNC__
2954 #define __FUNC__ /*<a name=""></a>*/"MatGetLocalSize"
2955 /*@
2956    MatGetLocalSize - Returns the number of rows and columns in a matrix
2957    stored locally.  This information may be implementation dependent, so
2958    use with care.
2959 
2960    Not Collective
2961 
2962    Input Parameters:
2963 .  mat - the matrix
2964 
2965    Output Parameters:
2966 +  m - the number of local rows
2967 -  n - the number of local columns
2968 
2969    Level: beginner
2970 
2971 .keywords: matrix, dimension, size, local, rows, columns, get
2972 
2973 .seealso: MatGetSize()
2974 @*/
2975 int MatGetLocalSize(Mat mat,int *m,int* n)
2976 {
2977   int ierr;
2978 
2979   PetscFunctionBegin;
2980   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2981   ierr = (*mat->ops->getlocalsize)(mat,m,n);CHKERRQ(ierr);
2982   PetscFunctionReturn(0);
2983 }
2984 
2985 #undef __FUNC__
2986 #define __FUNC__ /*<a name=""></a>*/"MatGetOwnershipRange"
2987 /*@
2988    MatGetOwnershipRange - Returns the range of matrix rows owned by
2989    this processor, assuming that the matrix is laid out with the first
2990    n1 rows on the first processor, the next n2 rows on the second, etc.
2991    For certain parallel layouts this range may not be well defined.
2992 
2993    Not Collective
2994 
2995    Input Parameters:
2996 .  mat - the matrix
2997 
2998    Output Parameters:
2999 +  m - the global index of the first local row
3000 -  n - one more than the global index of the last local row
3001 
3002    Level: beginner
3003 
3004 .keywords: matrix, get, range, ownership
3005 @*/
3006 int MatGetOwnershipRange(Mat mat,int *m,int* n)
3007 {
3008   int ierr;
3009 
3010   PetscFunctionBegin;
3011   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3012   if (m) PetscValidIntPointer(m);
3013   if (n) PetscValidIntPointer(n);
3014   if (!mat->ops->getownershiprange) SETERRQ(PETSC_ERR_SUP,0,"");
3015   ierr = (*mat->ops->getownershiprange)(mat,m,n);CHKERRQ(ierr);
3016   PetscFunctionReturn(0);
3017 }
3018 
3019 #undef __FUNC__
3020 #define __FUNC__ /*<a name=""></a>*/"MatILUFactorSymbolic"
3021 /*@
3022    MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
3023    Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric()
3024    to complete the factorization.
3025 
3026    Collective on Mat
3027 
3028    Input Parameters:
3029 +  mat - the matrix
3030 .  row - row permutation
3031 .  column - column permutation
3032 -  info - structure containing
3033 $      levels - number of levels of fill.
3034 $      expected fill - as ratio of original fill.
3035 $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
3036                 missing diagonal entries)
3037 
3038    Output Parameters:
3039 .  fact - new matrix that has been symbolically factored
3040 
3041    Notes:
3042    See the users manual for additional information about
3043    choosing the fill factor for better efficiency.
3044 
3045    Most users should employ the simplified SLES interface for linear solvers
3046    instead of working directly with matrix algebra routines such as this.
3047    See, e.g., SLESCreate().
3048 
3049    Level: developer
3050 
3051 .keywords: matrix, factor, incomplete, ILU, symbolic, fill
3052 
3053 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
3054           MatGetOrdering()
3055 
3056 @*/
3057 int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatILUInfo *info,Mat *fact)
3058 {
3059   int ierr;
3060 
3061   PetscFunctionBegin;
3062   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3063   PetscValidPointer(fact);
3064   if (info && info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Levels of fill negative %d",info->levels);
3065   if (info && info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",info->fill);
3066   if (!mat->ops->ilufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Only MatCreateMPIRowbs() matrices support parallel ILU");
3067   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3068   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3069 
3070   PLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);
3071   ierr = (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr);
3072   PLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);
3073   PetscFunctionReturn(0);
3074 }
3075 
3076 #undef __FUNC__
3077 #define __FUNC__ /*<a name=""></a>*/"MatIncompleteCholeskyFactorSymbolic"
3078 /*@
3079    MatIncompleteCholeskyFactorSymbolic - Performs symbolic incomplete
3080    Cholesky factorization for a symmetric matrix.  Use
3081    MatCholeskyFactorNumeric() to complete the factorization.
3082 
3083    Collective on Mat
3084 
3085    Input Parameters:
3086 +  mat - the matrix
3087 .  perm - row and column permutation
3088 .  fill - levels of fill
3089 -  f - expected fill as ratio of original fill
3090 
3091    Output Parameter:
3092 .  fact - the factored matrix
3093 
3094    Notes:
3095    Currently only no-fill factorization is supported.
3096 
3097    Most users should employ the simplified SLES interface for linear solvers
3098    instead of working directly with matrix algebra routines such as this.
3099    See, e.g., SLESCreate().
3100 
3101    Level: developer
3102 
3103 .keywords: matrix, factor, incomplete, ICC, Cholesky, symbolic, fill
3104 
3105 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor()
3106 @*/
3107 int MatIncompleteCholeskyFactorSymbolic(Mat mat,IS perm,PetscReal f,int fill,Mat *fact)
3108 {
3109   int ierr;
3110 
3111   PetscFunctionBegin;
3112   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3113   PetscValidPointer(fact);
3114   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3115   if (fill < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Fill negative %d",fill);
3116   if (f < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",f);
3117   if (!mat->ops->incompletecholeskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support ICC in parallel");
3118   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3119 
3120   PLogEventBegin(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);
3121   ierr = (*mat->ops->incompletecholeskyfactorsymbolic)(mat,perm,f,fill,fact);CHKERRQ(ierr);
3122   PLogEventEnd(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);
3123   PetscFunctionReturn(0);
3124 }
3125 
3126 #undef __FUNC__
3127 #define __FUNC__ /*<a name=""></a>*/"MatGetArray"
3128 /*@C
3129    MatGetArray - Returns a pointer to the element values in the matrix.
3130    The result of this routine is dependent on the underlying matrix data
3131    structure, and may not even work for certain matrix types.  You MUST
3132    call MatRestoreArray() when you no longer need to access the array.
3133 
3134    Not Collective
3135 
3136    Input Parameter:
3137 .  mat - the matrix
3138 
3139    Output Parameter:
3140 .  v - the location of the values
3141 
3142    Currently returns an array only for the dense formats, giving access to
3143    the local portion of the matrix in the usual Fortran column-oriented format.
3144 
3145    Fortran Note:
3146    This routine is used differently from Fortran, e.g.,
3147 .vb
3148         Mat         mat
3149         Scalar      mat_array(1)
3150         PetscOffset i_mat
3151         int         ierr
3152         call MatGetArray(mat,mat_array,i_mat,ierr)
3153 
3154   C  Access first local entry in matrix; note that array is
3155   C  treated as one dimensional
3156         value = mat_array(i_mat + 1)
3157 
3158         [... other code ...]
3159         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3160 .ve
3161 
3162    See the Fortran chapter of the users manual and
3163    petsc/src/mat/examples/tests for details.
3164 
3165    Level: advanced
3166 
3167 .keywords: matrix, array, elements, values
3168 
3169 .seealso: MatRestoreArray(), MatGetArrayF90()
3170 @*/
3171 int MatGetArray(Mat mat,Scalar **v)
3172 {
3173   int ierr;
3174 
3175   PetscFunctionBegin;
3176   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3177   PetscValidPointer(v);
3178   if (!mat->ops->getarray) SETERRQ(PETSC_ERR_SUP,0,"");
3179   ierr = (*mat->ops->getarray)(mat,v);CHKERRQ(ierr);
3180   PetscFunctionReturn(0);
3181 }
3182 
3183 #undef __FUNC__
3184 #define __FUNC__ /*<a name=""></a>*/"MatRestoreArray"
3185 /*@C
3186    MatRestoreArray - Restores the matrix after MatGetArray() has been called.
3187 
3188    Not Collective
3189 
3190    Input Parameter:
3191 +  mat - the matrix
3192 -  v - the location of the values
3193 
3194    Fortran Note:
3195    This routine is used differently from Fortran, e.g.,
3196 .vb
3197         Mat         mat
3198         Scalar      mat_array(1)
3199         PetscOffset i_mat
3200         int         ierr
3201         call MatGetArray(mat,mat_array,i_mat,ierr)
3202 
3203   C  Access first local entry in matrix; note that array is
3204   C  treated as one dimensional
3205         value = mat_array(i_mat + 1)
3206 
3207         [... other code ...]
3208         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3209 .ve
3210 
3211    See the Fortran chapter of the users manual and
3212    petsc/src/mat/examples/tests for details
3213 
3214    Level: advanced
3215 
3216 .keywords: matrix, array, elements, values, restore
3217 
3218 .seealso: MatGetArray(), MatRestoreArrayF90()
3219 @*/
3220 int MatRestoreArray(Mat mat,Scalar **v)
3221 {
3222   int ierr;
3223 
3224   PetscFunctionBegin;
3225   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3226   PetscValidPointer(v);
3227   if (!mat->ops->restorearray) SETERRQ(PETSC_ERR_SUP,0,"");
3228   ierr = (*mat->ops->restorearray)(mat,v);CHKERRQ(ierr);
3229   PetscFunctionReturn(0);
3230 }
3231 
3232 #undef __FUNC__
3233 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrices"
3234 /*@C
3235    MatGetSubMatrices - Extracts several submatrices from a matrix. If submat
3236    points to an array of valid matrices, they may be reused to store the new
3237    submatrices.
3238 
3239    Collective on Mat
3240 
3241    Input Parameters:
3242 +  mat - the matrix
3243 .  n   - the number of submatrixes to be extracted (on this processor, may be zero)
3244 .  irow, icol - index sets of rows and columns to extract
3245 -  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3246 
3247    Output Parameter:
3248 .  submat - the array of submatrices
3249 
3250    Notes:
3251    MatGetSubMatrices() can extract only sequential submatrices
3252    (from both sequential and parallel matrices). Use MatGetSubMatrix()
3253    to extract a parallel submatrix.
3254 
3255    When extracting submatrices from a parallel matrix, each processor can
3256    form a different submatrix by setting the rows and columns of its
3257    individual index sets according to the local submatrix desired.
3258 
3259    When finished using the submatrices, the user should destroy
3260    them with MatDestroySubMatrices().
3261 
3262    MAT_REUSE_MATRIX can only be used when the nonzero structure of the
3263    original matrix has not changed from that last call to MatGetSubMatrices().
3264 
3265    Fortran Note:
3266    The Fortran interface is slightly different from that given below; it
3267    requires one to pass in  as submat a Mat (integer) array of size at least m.
3268 
3269    Level: advanced
3270 
3271 .keywords: matrix, get, submatrix, submatrices
3272 
3273 .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal()
3274 @*/
3275 int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatReuse scall,Mat **submat)
3276 {
3277   int ierr;
3278 
3279   PetscFunctionBegin;
3280   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3281   if (!mat->ops->getsubmatrices) SETERRQ(PETSC_ERR_SUP,0,"");
3282   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3283 
3284   PLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);
3285   ierr = (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);CHKERRQ(ierr);
3286   PLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);
3287 
3288   PetscFunctionReturn(0);
3289 }
3290 
3291 #undef __FUNC__
3292 #define __FUNC__ /*<a name=""></a>*/"MatDestroyMatrices"
3293 /*@C
3294    MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices().
3295 
3296    Collective on Mat
3297 
3298    Input Parameters:
3299 +  n - the number of local matrices
3300 -  mat - the matrices
3301 
3302    Level: advanced
3303 
3304 .keywords: matrix, destroy, submatrix, submatrices
3305 
3306 .seealso: MatGetSubMatrices()
3307 @*/
3308 int MatDestroyMatrices(int n,Mat **mat)
3309 {
3310   int ierr,i;
3311 
3312   PetscFunctionBegin;
3313   if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,1,"Trying to destroy negative number of matrices %d",n);
3314   PetscValidPointer(mat);
3315   for (i=0; i<n; i++) {
3316     ierr = MatDestroy((*mat)[i]);CHKERRQ(ierr);
3317   }
3318   /* memory is allocated even if n = 0 */
3319   ierr = PetscFree(*mat);CHKERRQ(ierr);
3320   PetscFunctionReturn(0);
3321 }
3322 
3323 #undef __FUNC__
3324 #define __FUNC__ /*<a name=""></a>*/"MatIncreaseOverlap"
3325 /*@
3326    MatIncreaseOverlap - Given a set of submatrices indicated by index sets,
3327    replaces the index sets by larger ones that represent submatrices with
3328    additional overlap.
3329 
3330    Collective on Mat
3331 
3332    Input Parameters:
3333 +  mat - the matrix
3334 .  n   - the number of index sets
3335 .  is  - the array of pointers to index sets
3336 -  ov  - the additional overlap requested
3337 
3338    Level: developer
3339 
3340 .keywords: matrix, overlap, Schwarz
3341 
3342 .seealso: MatGetSubMatrices()
3343 @*/
3344 int MatIncreaseOverlap(Mat mat,int n,IS *is,int ov)
3345 {
3346   int ierr;
3347 
3348   PetscFunctionBegin;
3349   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3350   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3351   if (mat->factor)     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3352 
3353   if (!ov) PetscFunctionReturn(0);
3354   if (!mat->ops->increaseoverlap) SETERRQ(PETSC_ERR_SUP,0,"");
3355   PLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);
3356   ierr = (*mat->ops->increaseoverlap)(mat,n,is,ov);CHKERRQ(ierr);
3357   PLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);
3358   PetscFunctionReturn(0);
3359 }
3360 
3361 #undef __FUNC__
3362 #define __FUNC__ /*<a name=""></a>*/"MatPrintHelp"
3363 /*@
3364    MatPrintHelp - Prints all the options for the matrix.
3365 
3366    Collective on Mat
3367 
3368    Input Parameter:
3369 .  mat - the matrix
3370 
3371    Options Database Keys:
3372 +  -help - Prints matrix options
3373 -  -h - Prints matrix options
3374 
3375    Level: developer
3376 
3377 .keywords: mat, help
3378 
3379 .seealso: MatCreate(), MatCreateXXX()
3380 @*/
3381 int MatPrintHelp(Mat mat)
3382 {
3383   static PetscTruth called = PETSC_FALSE;
3384   int               ierr;
3385   MPI_Comm          comm;
3386 
3387   PetscFunctionBegin;
3388   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3389 
3390   comm = mat->comm;
3391   if (!called) {
3392     ierr = (*PetscHelpPrintf)(comm,"General matrix options:\n");CHKERRQ(ierr);
3393     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_info: view basic matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr);
3394     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_info_detailed: view detailed matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr);
3395     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_draw: draw nonzero matrix structure during MatAssemblyEnd()\n");CHKERRQ(ierr);
3396     ierr = (*PetscHelpPrintf)(comm,"      -draw_pause <sec>: set seconds of display pause\n");CHKERRQ(ierr);
3397     ierr = (*PetscHelpPrintf)(comm,"      -display <name>: set alternate display\n");CHKERRQ(ierr);
3398     called = PETSC_TRUE;
3399   }
3400   if (mat->ops->printhelp) {
3401     ierr = (*mat->ops->printhelp)(mat);CHKERRQ(ierr);
3402   }
3403   PetscFunctionReturn(0);
3404 }
3405 
3406 #undef __FUNC__
3407 #define __FUNC__ /*<a name=""></a>*/"MatGetBlockSize"
3408 /*@
3409    MatGetBlockSize - Returns the matrix block size; useful especially for the
3410    block row and block diagonal formats.
3411 
3412    Not Collective
3413 
3414    Input Parameter:
3415 .  mat - the matrix
3416 
3417    Output Parameter:
3418 .  bs - block size
3419 
3420    Notes:
3421    Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG.
3422    Block row formats are MATSEQBAIJ, MATMPIBAIJ
3423 
3424    Level: intermediate
3425 
3426 .keywords: matrix, get, block, size
3427 
3428 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag()
3429 @*/
3430 int MatGetBlockSize(Mat mat,int *bs)
3431 {
3432   int ierr;
3433 
3434   PetscFunctionBegin;
3435   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3436   PetscValidIntPointer(bs);
3437   if (!mat->ops->getblocksize) SETERRQ(PETSC_ERR_SUP,0,"");
3438   ierr = (*mat->ops->getblocksize)(mat,bs);CHKERRQ(ierr);
3439   PetscFunctionReturn(0);
3440 }
3441 
3442 #undef __FUNC__
3443 #define __FUNC__ /*<a name=""></a>*/"MatGetRowIJ"
3444 /*@C
3445     MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices.
3446 
3447    Collective on Mat
3448 
3449     Input Parameters:
3450 +   mat - the matrix
3451 .   shift -  0 or 1 indicating we want the indices starting at 0 or 1
3452 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3453                 symmetrized
3454 
3455     Output Parameters:
3456 +   n - number of rows in the (possibly compressed) matrix
3457 .   ia - the row pointers
3458 .   ja - the column indices
3459 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3460 
3461     Level: developer
3462 
3463 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3464 @*/
3465 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3466 {
3467   int ierr;
3468 
3469   PetscFunctionBegin;
3470   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3471   if (ia) PetscValidIntPointer(ia);
3472   if (ja) PetscValidIntPointer(ja);
3473   PetscValidIntPointer(done);
3474   if (!mat->ops->getrowij) *done = PETSC_FALSE;
3475   else {
3476     *done = PETSC_TRUE;
3477     ierr  = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3478   }
3479   PetscFunctionReturn(0);
3480 }
3481 
3482 #undef __FUNC__
3483 #define __FUNC__ /*<a name=""></a>*/"MatGetColumnIJ"
3484 /*@C
3485     MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices.
3486 
3487     Collective on Mat
3488 
3489     Input Parameters:
3490 +   mat - the matrix
3491 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3492 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3493                 symmetrized
3494 
3495     Output Parameters:
3496 +   n - number of columns in the (possibly compressed) matrix
3497 .   ia - the column pointers
3498 .   ja - the row indices
3499 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3500 
3501     Level: developer
3502 
3503 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3504 @*/
3505 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3506 {
3507   int ierr;
3508 
3509   PetscFunctionBegin;
3510   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3511   if (ia) PetscValidIntPointer(ia);
3512   if (ja) PetscValidIntPointer(ja);
3513   PetscValidIntPointer(done);
3514 
3515   if (!mat->ops->getcolumnij) *done = PETSC_FALSE;
3516   else {
3517     *done = PETSC_TRUE;
3518     ierr  = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3519   }
3520   PetscFunctionReturn(0);
3521 }
3522 
3523 #undef __FUNC__
3524 #define __FUNC__ /*<a name=""></a>*/"MatRestoreRowIJ"
3525 /*@C
3526     MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with
3527     MatGetRowIJ().
3528 
3529     Collective on Mat
3530 
3531     Input Parameters:
3532 +   mat - the matrix
3533 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3534 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3535                 symmetrized
3536 
3537     Output Parameters:
3538 +   n - size of (possibly compressed) matrix
3539 .   ia - the row pointers
3540 .   ja - the column indices
3541 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3542 
3543     Level: developer
3544 
3545 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3546 @*/
3547 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3548 {
3549   int ierr;
3550 
3551   PetscFunctionBegin;
3552   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3553   if (ia) PetscValidIntPointer(ia);
3554   if (ja) PetscValidIntPointer(ja);
3555   PetscValidIntPointer(done);
3556 
3557   if (!mat->ops->restorerowij) *done = PETSC_FALSE;
3558   else {
3559     *done = PETSC_TRUE;
3560     ierr  = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3561   }
3562   PetscFunctionReturn(0);
3563 }
3564 
3565 #undef __FUNC__
3566 #define __FUNC__ /*<a name=""></a>*/"MatRestoreColumnIJ"
3567 /*@C
3568     MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with
3569     MatGetColumnIJ().
3570 
3571     Collective on Mat
3572 
3573     Input Parameters:
3574 +   mat - the matrix
3575 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3576 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3577                 symmetrized
3578 
3579     Output Parameters:
3580 +   n - size of (possibly compressed) matrix
3581 .   ia - the column pointers
3582 .   ja - the row indices
3583 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3584 
3585     Level: developer
3586 
3587 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3588 @*/
3589 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3590 {
3591   int ierr;
3592 
3593   PetscFunctionBegin;
3594   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3595   if (ia) PetscValidIntPointer(ia);
3596   if (ja) PetscValidIntPointer(ja);
3597   PetscValidIntPointer(done);
3598 
3599   if (!mat->ops->restorecolumnij) *done = PETSC_FALSE;
3600   else {
3601     *done = PETSC_TRUE;
3602     ierr  = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3603   }
3604   PetscFunctionReturn(0);
3605 }
3606 
3607 #undef __FUNC__
3608 #define __FUNC__ /*<a name=""></a>*/"MatColoringPatch"
3609 /*@C
3610     MatColoringPatch -Used inside matrix coloring routines that
3611     use MatGetRowIJ() and/or MatGetColumnIJ().
3612 
3613     Collective on Mat
3614 
3615     Input Parameters:
3616 +   mat - the matrix
3617 .   n   - number of colors
3618 -   colorarray - array indicating color for each column
3619 
3620     Output Parameters:
3621 .   iscoloring - coloring generated using colorarray information
3622 
3623     Level: developer
3624 
3625 .seealso: MatGetRowIJ(), MatGetColumnIJ()
3626 
3627 .keywords: mat, coloring, patch
3628 @*/
3629 int MatColoringPatch(Mat mat,int n,int *colorarray,ISColoring *iscoloring)
3630 {
3631   int ierr;
3632 
3633   PetscFunctionBegin;
3634   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3635   PetscValidIntPointer(colorarray);
3636 
3637   if (!mat->ops->coloringpatch) {SETERRQ(PETSC_ERR_SUP,0,"");}
3638   else {
3639     ierr  = (*mat->ops->coloringpatch)(mat,n,colorarray,iscoloring);CHKERRQ(ierr);
3640   }
3641   PetscFunctionReturn(0);
3642 }
3643 
3644 
3645 #undef __FUNC__
3646 #define __FUNC__ /*<a name=""></a>*/"MatSetUnfactored"
3647 /*@
3648    MatSetUnfactored - Resets a factored matrix to be treated as unfactored.
3649 
3650    Collective on Mat
3651 
3652    Input Parameter:
3653 .  mat - the factored matrix to be reset
3654 
3655    Notes:
3656    This routine should be used only with factored matrices formed by in-place
3657    factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE
3658    format).  This option can save memory, for example, when solving nonlinear
3659    systems with a matrix-free Newton-Krylov method and a matrix-based, in-place
3660    ILU(0) preconditioner.
3661 
3662    Note that one can specify in-place ILU(0) factorization by calling
3663 .vb
3664      PCType(pc,PCILU);
3665      PCILUSeUseInPlace(pc);
3666 .ve
3667    or by using the options -pc_type ilu -pc_ilu_in_place
3668 
3669    In-place factorization ILU(0) can also be used as a local
3670    solver for the blocks within the block Jacobi or additive Schwarz
3671    methods (runtime option: -sub_pc_ilu_in_place).  See the discussion
3672    of these preconditioners in the users manual for details on setting
3673    local solver options.
3674 
3675    Most users should employ the simplified SLES interface for linear solvers
3676    instead of working directly with matrix algebra routines such as this.
3677    See, e.g., SLESCreate().
3678 
3679    Level: developer
3680 
3681 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace()
3682 
3683 .keywords: matrix-free, in-place ILU, in-place LU
3684 @*/
3685 int MatSetUnfactored(Mat mat)
3686 {
3687   int ierr;
3688 
3689   PetscFunctionBegin;
3690   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3691   mat->factor = 0;
3692   if (!mat->ops->setunfactored) PetscFunctionReturn(0);
3693   ierr = (*mat->ops->setunfactored)(mat);CHKERRQ(ierr);
3694   PetscFunctionReturn(0);
3695 }
3696 
3697 #undef __FUNC__
3698 #define __FUNC__ /*<a name=""></a>*/"MatGetType"
3699 /*@C
3700    MatGetType - Gets the matrix type and name (as a string) from the matrix.
3701 
3702    Not Collective
3703 
3704    Input Parameter:
3705 .  mat - the matrix
3706 
3707    Output Parameter:
3708 +  type - the matrix type (or use PETSC_NULL)
3709 -  name - name of matrix type (or use PETSC_NULL)
3710 
3711    Level: intermediate
3712 
3713 .keywords: matrix, get, type, name
3714 @*/
3715 int MatGetType(Mat mat,MatType *type,char **name)
3716 {
3717   int  itype = (int)mat->type;
3718   char *matname[10];
3719 
3720   PetscFunctionBegin;
3721   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3722 
3723   if (type) *type = (MatType) mat->type;
3724   if (name) {
3725     /* Note:  Be sure that this list corresponds to the enum in petscmat.h */
3726     matname[0] = "MATSEQDENSE";
3727     matname[1] = "MATSEQAIJ";
3728     matname[2] = "MATMPIAIJ";
3729     matname[3] = "MATSHELL";
3730     matname[4] = "MATMPIROWBS";
3731     matname[5] = "MATSEQBDIAG";
3732     matname[6] = "MATMPIBDIAG";
3733     matname[7] = "MATMPIDENSE";
3734     matname[8] = "MATSEQBAIJ";
3735     matname[9] = "MATMPIBAIJ";
3736 
3737     if (itype < 0 || itype > 9) *name = "Unknown matrix type";
3738     else                        *name = matname[itype];
3739   }
3740   PetscFunctionReturn(0);
3741 }
3742 
3743 /*MC
3744     MatGetArrayF90 - Accesses a matrix array from Fortran90.
3745 
3746     Synopsis:
3747     MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3748 
3749     Not collective
3750 
3751     Input Parameter:
3752 .   x - matrix
3753 
3754     Output Parameters:
3755 +   xx_v - the Fortran90 pointer to the array
3756 -   ierr - error code
3757 
3758     Example of Usage:
3759 .vb
3760       Scalar, pointer xx_v(:)
3761       ....
3762       call MatGetArrayF90(x,xx_v,ierr)
3763       a = xx_v(3)
3764       call MatRestoreArrayF90(x,xx_v,ierr)
3765 .ve
3766 
3767     Notes:
3768     Not yet supported for all F90 compilers
3769 
3770     Level: advanced
3771 
3772 .seealso:  MatRestoreArrayF90(), MatGetArray(), MatRestoreArray()
3773 
3774 .keywords:  matrix, array, f90
3775 M*/
3776 
3777 /*MC
3778     MatRestoreArrayF90 - Restores a matrix array that has been
3779     accessed with MatGetArrayF90().
3780 
3781     Synopsis:
3782     MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3783 
3784     Not collective
3785 
3786     Input Parameters:
3787 +   x - matrix
3788 -   xx_v - the Fortran90 pointer to the array
3789 
3790     Output Parameter:
3791 .   ierr - error code
3792 
3793     Example of Usage:
3794 .vb
3795        Scalar, pointer xx_v(:)
3796        ....
3797        call MatGetArrayF90(x,xx_v,ierr)
3798        a = xx_v(3)
3799        call MatRestoreArrayF90(x,xx_v,ierr)
3800 .ve
3801 
3802     Notes:
3803     Not yet supported for all F90 compilers
3804 
3805     Level: advanced
3806 
3807 .seealso:  MatGetArrayF90(), MatGetArray(), MatRestoreArray()
3808 
3809 .keywords:  matrix, array, f90
3810 M*/
3811 
3812 
3813 #undef __FUNC__
3814 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrix"
3815 /*@
3816     MatGetSubMatrix - Gets a single submatrix on the same number of processors
3817                       as the original matrix.
3818 
3819     Collective on Mat
3820 
3821     Input Parameters:
3822 +   mat - the original matrix
3823 .   isrow - rows this processor should obtain
3824 .   iscol - columns for all processors you wish to keep
3825 .   csize - number of columns "local" to this processor (does nothing for sequential
3826             matrices). This should match the result from VecGetLocalSize(x,...) if you
3827             plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE
3828 -   cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3829 
3830     Output Parameter:
3831 .   newmat - the new submatrix, of the same type as the old
3832 
3833     Level: advanced
3834 
3835 .keywords: matrix, get, submatrix, submatrices
3836 
3837 .seealso: MatGetSubMatrices()
3838 @*/
3839 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat)
3840 {
3841   int     ierr, size;
3842   Mat     *local;
3843 
3844   PetscFunctionBegin;
3845   ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr);
3846 
3847   /* if original matrix is on just one processor then use submatrix generated */
3848   if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) {
3849     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr);
3850     PetscFunctionReturn(0);
3851   } else if (!mat->ops->getsubmatrix && size == 1) {
3852     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr);
3853     *newmat = *local;
3854     ierr = PetscFree(local);CHKERRQ(ierr);
3855     PetscFunctionReturn(0);
3856   }
3857 
3858   if (!mat->ops->getsubmatrix) SETERRQ(PETSC_ERR_SUP,0,"Not currently implemented");
3859   ierr = (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr);
3860   PetscFunctionReturn(0);
3861 }
3862 
3863 #undef __FUNC__
3864 #define __FUNC__ /*<a name=""></a>*/"MatGetMaps"
3865 /*@C
3866    MatGetMaps - Returns the maps associated with the matrix.
3867 
3868    Not Collective
3869 
3870    Input Parameter:
3871 .  mat - the matrix
3872 
3873    Output Parameters:
3874 +  rmap - the row (right) map
3875 -  cmap - the column (left) map
3876 
3877    Level: developer
3878 
3879 .keywords: matrix, get, map
3880 @*/
3881 int MatGetMaps(Mat mat,Map *rmap,Map *cmap)
3882 {
3883   int ierr;
3884 
3885   PetscFunctionBegin;
3886   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3887   ierr = (*mat->ops->getmaps)(mat,rmap,cmap);CHKERRQ(ierr);
3888   PetscFunctionReturn(0);
3889 }
3890 
3891 /*
3892       Version that works for all PETSc matrices
3893 */
3894 #undef __FUNC__
3895 #define __FUNC__ /*<a name=""></a>*/"MatGetMaps_Petsc"
3896 int MatGetMaps_Petsc(Mat mat,Map *rmap,Map *cmap)
3897 {
3898   PetscFunctionBegin;
3899   if (rmap) *rmap = mat->rmap;
3900   if (cmap) *cmap = mat->cmap;
3901   PetscFunctionReturn(0);
3902 }
3903 
3904 #undef __FUNC__
3905 #define __FUNC__ /*<a name=""></a>*/"MatSetStashInitialSize"
3906 /*@
3907    MatSetStashInitialSize - sets the sizes of the matrix stash, that is
3908    used during the assembly process to store values that belong to
3909    other processors.
3910 
3911    Not Collective
3912 
3913    Input Parameters:
3914 +  mat   - the matrix
3915 .  size  - the initial size of the stash.
3916 -  bsize - the initial size of the block-stash(if used).
3917 
3918    Options Database Keys:
3919 +   -matstash_initial_size <size> or <size0,size1,...sizep-1>
3920 -   -matstash_block_initial_size <bsize>  or <bsize0,bsize1,...bsizep-1>
3921 
3922    Level: intermediate
3923 
3924    Notes:
3925      The block-stash is used for values set with VecSetValuesBlocked() while
3926      the stash is used for values set with VecSetValues()
3927 
3928      Run with the option -log_info and look for output of the form
3929      MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
3930      to determine the appropriate value, MM, to use for size and
3931      MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
3932      to determine the value, BMM to use for bsize
3933 
3934 .keywords: matrix, stash, assembly
3935 @*/
3936 int MatSetStashInitialSize(Mat mat,int size, int bsize)
3937 {
3938   int ierr;
3939 
3940   PetscFunctionBegin;
3941   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3942   ierr = MatStashSetInitialSize_Private(&mat->stash,size);CHKERRQ(ierr);
3943   ierr = MatStashSetInitialSize_Private(&mat->bstash,bsize);CHKERRQ(ierr);
3944   PetscFunctionReturn(0);
3945 }
3946 
3947 #undef __FUNC__
3948 #define __FUNC__ /*<a name=""></a>*/"MatInterpolateAdd"
3949 /*@
3950    MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of
3951      the matrix
3952 
3953    Collective on Mat
3954 
3955    Input Parameters:
3956 +  mat   - the matrix
3957 .  x,y - the vectors
3958 -  w - where the result is stored
3959 
3960    Level: intermediate
3961 
3962    Notes:
3963     w may be the same vector as y.
3964 
3965     This allows one to use either the restriction or interpolation (its transpose)
3966     matrix to do the interpolation
3967 
3968 .keywords: interpolate,
3969 
3970 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
3971 
3972 @*/
3973 int MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w)
3974 {
3975   int M,N,ierr;
3976 
3977   PetscFunctionBegin;
3978   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
3979   if (N > M) {
3980     ierr = MatMultTransposeAdd(A,x,y,w);CHKERRQ(ierr);
3981   } else {
3982     ierr = MatMultAdd(A,x,y,w);CHKERRQ(ierr);
3983   }
3984   PetscFunctionReturn(0);
3985 }
3986 
3987 #undef __FUNC__
3988 #define __FUNC__ /*<a name=""></a>*/"MatInterpolate"
3989 /*@
3990    MatInterpolate - y = A*x or A'*x depending on the shape of
3991      the matrix
3992 
3993    Collective on Mat
3994 
3995    Input Parameters:
3996 +  mat   - the matrix
3997 -  x,y - the vectors
3998 
3999    Level: intermediate
4000 
4001    Notes:
4002     This allows one to use either the restriction or interpolation (its transpose)
4003     matrix to do the interpolation
4004 
4005 .keywords: interpolate,
4006 
4007 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
4008 
4009 @*/
4010 int MatInterpolate(Mat A,Vec x,Vec y)
4011 {
4012   int M,N,ierr;
4013 
4014   PetscFunctionBegin;
4015   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
4016   if (N > M) {
4017     ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
4018   } else {
4019     ierr = MatMult(A,x,y);CHKERRQ(ierr);
4020   }
4021   PetscFunctionReturn(0);
4022 }
4023 
4024 #undef __FUNC__
4025 #define __FUNC__ /*<a name=""></a>*/"MatRestrict"
4026 /*@
4027    MatRestrict - y = A*x or A'*x
4028 
4029    Collective on Mat
4030 
4031    Input Parameters:
4032 +  mat   - the matrix
4033 -  x,y - the vectors
4034 
4035    Level: intermediate
4036 
4037    Notes:
4038     This allows one to use either the restriction or interpolation (its transpose)
4039     matrix to do the restriction
4040 
4041 .keywords: interpolate,
4042 
4043 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatInterpolate()
4044 
4045 @*/
4046 int MatRestrict(Mat A,Vec x,Vec y)
4047 {
4048   int M,N,ierr;
4049 
4050   PetscFunctionBegin;
4051   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
4052   if (N > M) {
4053     ierr = MatMult(A,x,y);CHKERRQ(ierr);
4054   } else {
4055     ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
4056   }
4057   PetscFunctionReturn(0);
4058 }
4059