xref: /petsc/src/mat/interface/matrix.c (revision 273d9f13de75c4ed17021f7f2c11eebb99d26f0d)
1 /*$Id: matrix.c,v 1.382 2000/09/28 21:10:52 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    Concepts: matrices^row access
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   PetscValidType(mat);
80   MatPreallocated(mat);
81   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
82   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
83   if (!mat->ops->getrow) SETERRQ(PETSC_ERR_SUP,"");
84   ierr = PLogEventBegin(MAT_GetRow,mat,0,0,0);CHKERRQ(ierr);
85   ierr = (*mat->ops->getrow)(mat,row,ncols,cols,vals);CHKERRQ(ierr);
86   ierr = PLogEventEnd(MAT_GetRow,mat,0,0,0);CHKERRQ(ierr);
87   PetscFunctionReturn(0);
88 }
89 
90 #undef __FUNC__
91 #define __FUNC__ /*<a name=""></a>*/"MatRestoreRow"
92 /*@C
93    MatRestoreRow - Frees any temporary space allocated by MatGetRow().
94 
95    Not Collective
96 
97    Input Parameters:
98 +  mat - the matrix
99 .  row - the row to get
100 .  ncols, cols - the number of nonzeros and their columns
101 -  vals - if nonzero the column values
102 
103    Notes:
104    This routine should be called after you have finished examining the entries.
105 
106    Fortran Notes:
107    The calling sequence from Fortran is
108 .vb
109    MatRestoreRow(matrix,row,ncols,cols,values,ierr)
110       Mat     matrix (input)
111       integer row    (input)
112       integer ncols  (output)
113       integer cols(maxcols) (output)
114       double precision (or double complex) values(maxcols) output
115 .ve
116    Where maxcols >= maximum nonzeros in any row of the matrix.
117 
118    In Fortran MatRestoreRow() MUST be called after MatGetRow()
119    before another call to MatGetRow() can be made.
120 
121    Level: advanced
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,"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    Concepts: matrices^viewing
186    Concepts: matrices^plotting
187    Concepts: matrices^printing
188 
189 .seealso: ViewerSetFormat(), ViewerASCIIOpen(), ViewerDrawOpen(),
190           ViewerSocketOpen(), ViewerBinaryOpen(), MatLoad()
191 @*/
192 int MatView(Mat mat,Viewer viewer)
193 {
194   int        format,ierr,rows,cols;
195   PetscTruth isascii;
196   char       *cstr;
197 
198   PetscFunctionBegin;
199   PetscValidHeaderSpecific(mat,MAT_COOKIE);
200   PetscValidType(mat);
201   MatPreallocated(mat);
202   if (!viewer) viewer = VIEWER_STDOUT_(mat->comm);
203   PetscValidHeaderSpecific(viewer,VIEWER_COOKIE);
204   PetscCheckSameComm(mat,viewer);
205   if (!mat->assembled) SETERRQ(1,"Must call MatAssemblyBegin/End() before viewing matrix");
206 
207   ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr);
208   if (isascii) {
209     ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr);
210     if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) {
211       ierr = ViewerASCIIPrintf(viewer,"Matrix Object:\n");CHKERRQ(ierr);
212       ierr = ViewerASCIIPushTab(viewer);CHKERRQ(ierr);
213       ierr = MatGetType(mat,&cstr);CHKERRQ(ierr);
214       ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr);
215       ierr = ViewerASCIIPrintf(viewer,"type=%s, rows=%d, cols=%d\n",cstr,rows,cols);CHKERRQ(ierr);
216       if (mat->ops->getinfo) {
217         MatInfo info;
218         ierr = MatGetInfo(mat,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr);
219         ierr = ViewerASCIIPrintf(viewer,"total: nonzeros=%d, allocated nonzeros=%d\n",
220                           (int)info.nz_used,(int)info.nz_allocated);CHKERRQ(ierr);
221       }
222     }
223   }
224   if (mat->ops->view) {
225     ierr = ViewerASCIIPushTab(viewer);CHKERRQ(ierr);
226     ierr = (*mat->ops->view)(mat,viewer);CHKERRQ(ierr);
227     ierr = ViewerASCIIPopTab(viewer);CHKERRQ(ierr);
228   } else if (!isascii) {
229     SETERRQ1(1,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
230   }
231   if (isascii) {
232     ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr);
233     if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) {
234       ierr = ViewerASCIIPopTab(viewer);CHKERRQ(ierr);
235     }
236   }
237   PetscFunctionReturn(0);
238 }
239 
240 #undef __FUNC__
241 #define __FUNC__ /*<a name=""></a>*/"MatScaleSystem"
242 /*@C
243    MatScaleSystem - Scale a vector solution and right hand side to
244    match the scaling of a scaled matrix.
245 
246    Collective on Mat
247 
248    Input Parameter:
249 +  mat - the matrix
250 .  x - solution vector (or PETSC_NULL)
251 +  b - right hand side vector (or PETSC_NULL)
252 
253 
254    Notes:
255    For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
256    internally scaled, so this does nothing. For MPIROWBS it
257    permutes and diagonally scales.
258 
259    The SLES methods automatically call this routine when required
260    (via PCPreSolve()) so it is rarely used directly.
261 
262    Level: Developer
263 
264    Concepts: matrices^scaling
265 
266 .seealso: MatUseScaledForm(), MatUnScaleSystem()
267 @*/
268 int MatScaleSystem(Mat mat,Vec x,Vec b)
269 {
270   int ierr;
271 
272   PetscFunctionBegin;
273   PetscValidHeaderSpecific(mat,MAT_COOKIE);
274   PetscValidType(mat);
275   MatPreallocated(mat);
276   if (x) {PetscValidHeaderSpecific(x,VEC_COOKIE);PetscCheckSameComm(mat,x);}
277   if (b) {PetscValidHeaderSpecific(b,VEC_COOKIE);PetscCheckSameComm(mat,b);}
278 
279   if (mat->ops->scalesystem) {
280     ierr = (*mat->ops->scalesystem)(mat,x,b);CHKERRQ(ierr);
281   }
282   PetscFunctionReturn(0);
283 }
284 
285 #undef __FUNC__
286 #define __FUNC__ /*<a name=""></a>*/"MatUnScaleSystem"
287 /*@C
288    MatUnScaleSystem - Unscales a vector solution and right hand side to
289    match the original scaling of a scaled matrix.
290 
291    Collective on Mat
292 
293    Input Parameter:
294 +  mat - the matrix
295 .  x - solution vector (or PETSC_NULL)
296 +  b - right hand side vector (or PETSC_NULL)
297 
298 
299    Notes:
300    For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
301    internally scaled, so this does nothing. For MPIROWBS it
302    permutes and diagonally scales.
303 
304    The SLES methods automatically call this routine when required
305    (via PCPreSolve()) so it is rarely used directly.
306 
307    Level: Developer
308 
309 .seealso: MatUseScaledForm(), MatScaleSystem()
310 @*/
311 int MatUnScaleSystem(Mat mat,Vec x,Vec b)
312 {
313   int ierr;
314 
315   PetscFunctionBegin;
316   PetscValidHeaderSpecific(mat,MAT_COOKIE);
317   PetscValidType(mat);
318   MatPreallocated(mat);
319   if (x) {PetscValidHeaderSpecific(x,VEC_COOKIE);PetscCheckSameComm(mat,x);}
320   if (b) {PetscValidHeaderSpecific(b,VEC_COOKIE);PetscCheckSameComm(mat,b);}
321   if (mat->ops->unscalesystem) {
322     ierr = (*mat->ops->unscalesystem)(mat,x,b);CHKERRQ(ierr);
323   }
324   PetscFunctionReturn(0);
325 }
326 
327 #undef __FUNC__
328 #define __FUNC__ /*<a name=""></a>*/"MatUseScaledForm"
329 /*@C
330    MatUseScaledForm - For matrix storage formats that scale the
331    matrix (for example MPIRowBS matrices are diagonally scaled on
332    assembly) indicates matrix operations (MatMult() etc) are
333    applied using the scaled matrix.
334 
335    Collective on Mat
336 
337    Input Parameter:
338 +  mat - the matrix
339 -  scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for
340             applying the original matrix
341 
342    Notes:
343    For scaled matrix formats, applying the original, unscaled matrix
344    will be slightly more expensive
345 
346    Level: Developer
347 
348 .seealso: MatScaleSystem(), MatUnScaleSystem()
349 @*/
350 int MatUseScaledForm(Mat mat,PetscTruth scaled)
351 {
352   int ierr;
353 
354   PetscFunctionBegin;
355   PetscValidHeaderSpecific(mat,MAT_COOKIE);
356   PetscValidType(mat);
357   MatPreallocated(mat);
358   if (mat->ops->usescaledform) {
359     ierr = (*mat->ops->usescaledform)(mat,scaled);CHKERRQ(ierr);
360   }
361   PetscFunctionReturn(0);
362 }
363 
364 #undef __FUNC__
365 #define __FUNC__ /*<a name=""></a>*/"MatDestroy"
366 /*@C
367    MatDestroy - Frees space taken by a matrix.
368 
369    Collective on Mat
370 
371    Input Parameter:
372 .  A - the matrix
373 
374    Level: beginner
375 
376 @*/
377 int MatDestroy(Mat A)
378 {
379   int ierr;
380 
381   PetscFunctionBegin;
382   PetscValidHeaderSpecific(A,MAT_COOKIE);
383   PetscValidType(A);
384   MatPreallocated(A);
385   if (--A->refct > 0) PetscFunctionReturn(0);
386 
387   /* if memory was published with AMS then destroy it */
388   ierr = PetscObjectDepublish(A);CHKERRQ(ierr);
389   if (A->mapping) {
390     ierr = ISLocalToGlobalMappingDestroy(A->mapping);CHKERRQ(ierr);
391   }
392   if (A->bmapping) {
393     ierr = ISLocalToGlobalMappingDestroy(A->bmapping);CHKERRQ(ierr);
394   }
395   if (A->rmap) {
396     ierr = MapDestroy(A->rmap);CHKERRQ(ierr);
397   }
398   if (A->cmap) {
399     ierr = MapDestroy(A->cmap);CHKERRQ(ierr);
400   }
401 
402   ierr = (*A->ops->destroy)(A);CHKERRQ(ierr);
403   PLogObjectDestroy(A);
404   PetscHeaderDestroy(A);
405   PetscFunctionReturn(0);
406 }
407 
408 #undef __FUNC__
409 #define __FUNC__ /*<a name=""></a>*/"MatValid"
410 /*@
411    MatValid - Checks whether a matrix object is valid.
412 
413    Collective on Mat
414 
415    Input Parameter:
416 .  m - the matrix to check
417 
418    Output Parameter:
419    flg - flag indicating matrix status, either
420    PETSC_TRUE if matrix is valid, or PETSC_FALSE otherwise.
421 
422    Level: developer
423 
424    Concepts: matrices^validity
425 @*/
426 int MatValid(Mat m,PetscTruth *flg)
427 {
428   PetscFunctionBegin;
429   PetscValidIntPointer(flg);
430   if (!m)                           *flg = PETSC_FALSE;
431   else if (m->cookie != MAT_COOKIE) *flg = PETSC_FALSE;
432   else                              *flg = PETSC_TRUE;
433   PetscFunctionReturn(0);
434 }
435 
436 #undef __FUNC__
437 #define __FUNC__ /*<a name=""></a>*/"MatSetValues"
438 /*@
439    MatSetValues - Inserts or adds a block of values into a matrix.
440    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
441    MUST be called after all calls to MatSetValues() have been completed.
442 
443    Not Collective
444 
445    Input Parameters:
446 +  mat - the matrix
447 .  v - a logically two-dimensional array of values
448 .  m, idxm - the number of rows and their global indices
449 .  n, idxn - the number of columns and their global indices
450 -  addv - either ADD_VALUES or INSERT_VALUES, where
451    ADD_VALUES adds values to any existing entries, and
452    INSERT_VALUES replaces existing entries with new values
453 
454    Notes:
455    By default the values, v, are row-oriented and unsorted.
456    See MatSetOption() for other options.
457 
458    Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES
459    options cannot be mixed without intervening calls to the assembly
460    routines.
461 
462    MatSetValues() uses 0-based row and column numbers in Fortran
463    as well as in C.
464 
465    Negative indices may be passed in idxm and idxn, these rows and columns are
466    simply ignored. This allows easily inserting element stiffness matrices
467    with homogeneous Dirchlet boundary conditions that you don't want represented
468    in the matrix.
469 
470    Efficiency Alert:
471    The routine MatSetValuesBlocked() may offer much better efficiency
472    for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
473 
474    Level: beginner
475 
476    Concepts: matrices^putting entries in
477 
478 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
479 @*/
480 int MatSetValues(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v,InsertMode addv)
481 {
482   int ierr;
483 
484   PetscFunctionBegin;
485   if (!m || !n) PetscFunctionReturn(0); /* no values to insert */
486   PetscValidHeaderSpecific(mat,MAT_COOKIE);
487   PetscValidType(mat);
488   MatPreallocated(mat);
489   PetscValidIntPointer(idxm);
490   PetscValidIntPointer(idxn);
491   PetscValidScalarPointer(v);
492   if (mat->insertmode == NOT_SET_VALUES) {
493     mat->insertmode = addv;
494   }
495 #if defined(PETSC_USE_BOPT_g)
496   else if (mat->insertmode != addv) {
497     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
498   }
499   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
500 #endif
501 
502   if (mat->assembled) {
503     mat->was_assembled = PETSC_TRUE;
504     mat->assembled     = PETSC_FALSE;
505   }
506   ierr = PLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr);
507   if (!mat->ops->setvalues) SETERRQ(PETSC_ERR_SUP,"Not supported for this matrix type");
508   ierr = (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr);
509   ierr = PLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr);
510   PetscFunctionReturn(0);
511 }
512 
513 #undef __FUNC__
514 #define __FUNC__ /*<a name=""></a>*/"MatSetValuesBlocked"
515 /*@
516    MatSetValuesBlocked - Inserts or adds a block of values into a matrix.
517 
518    Not Collective
519 
520    Input Parameters:
521 +  mat - the matrix
522 .  v - a logically two-dimensional array of values
523 .  m, idxm - the number of block rows and their global block indices
524 .  n, idxn - the number of block columns and their global block indices
525 -  addv - either ADD_VALUES or INSERT_VALUES, where
526    ADD_VALUES adds values to any existing entries, and
527    INSERT_VALUES replaces existing entries with new values
528 
529    Notes:
530    By default the values, v, are row-oriented and unsorted. So the layout of
531    v is the same as for MatSetValues(). See MatSetOption() for other options.
532 
533    Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
534    options cannot be mixed without intervening calls to the assembly
535    routines.
536 
537    MatSetValuesBlocked() uses 0-based row and column numbers in Fortran
538    as well as in C.
539 
540    Negative indices may be passed in idxm and idxn, these rows and columns are
541    simply ignored. This allows easily inserting element stiffness matrices
542    with homogeneous Dirchlet boundary conditions that you don't want represented
543    in the matrix.
544 
545    Each time an entry is set within a sparse matrix via MatSetValues(),
546    internal searching must be done to determine where to place the the
547    data in the matrix storage space.  By instead inserting blocks of
548    entries via MatSetValuesBlocked(), the overhead of matrix assembly is
549    reduced.
550 
551    Restrictions:
552    MatSetValuesBlocked() is currently supported only for the block AIJ
553    matrix format (MATSEQBAIJ and MATMPIBAIJ, which are created via
554    MatCreateSeqBAIJ() and MatCreateMPIBAIJ()).
555 
556    Level: intermediate
557 
558    Concepts: matrices^putting entries in blocked
559 
560 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
561 @*/
562 int MatSetValuesBlocked(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v,InsertMode addv)
563 {
564   int ierr;
565 
566   PetscFunctionBegin;
567   if (!m || !n) PetscFunctionReturn(0); /* no values to insert */
568   PetscValidHeaderSpecific(mat,MAT_COOKIE);
569   PetscValidType(mat);
570   MatPreallocated(mat);
571   PetscValidIntPointer(idxm);
572   PetscValidIntPointer(idxn);
573   PetscValidScalarPointer(v);
574   if (mat->insertmode == NOT_SET_VALUES) {
575     mat->insertmode = addv;
576   }
577 #if defined(PETSC_USE_BOPT_g)
578   else if (mat->insertmode != addv) {
579     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
580   }
581   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
582 #endif
583 
584   if (mat->assembled) {
585     mat->was_assembled = PETSC_TRUE;
586     mat->assembled     = PETSC_FALSE;
587   }
588   ierr = PLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr);
589   if (!mat->ops->setvaluesblocked) SETERRQ(PETSC_ERR_SUP,"Not supported for this matrix type");
590   ierr = (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr);
591   ierr = PLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr);
592   PetscFunctionReturn(0);
593 }
594 
595 /*MC
596    MatSetValue - Set a single entry into a matrix.
597 
598    Synopsis:
599    void MatSetValue(Mat m,int row,int col,Scalar value,InsertMode mode);
600 
601    Not collective
602 
603    Input Parameters:
604 +  m - the matrix
605 .  row - the row location of the entry
606 .  col - the column location of the entry
607 .  value - the value to insert
608 -  mode - either INSERT_VALUES or ADD_VALUES
609 
610    Notes:
611    For efficiency one should use MatSetValues() and set several or many
612    values simultaneously if possible.
613 
614    Note that VecSetValue() does NOT return an error code (since this
615    is checked internally).
616 
617    Level: beginner
618 
619 .seealso: MatSetValues()
620 M*/
621 
622 #undef __FUNC__
623 #define __FUNC__ /*<a name=""></a>*/"MatGetValues"
624 /*@
625    MatGetValues - Gets a block of values from a matrix.
626 
627    Not Collective; currently only returns a local block
628 
629    Input Parameters:
630 +  mat - the matrix
631 .  v - a logically two-dimensional array for storing the values
632 .  m, idxm - the number of rows and their global indices
633 -  n, idxn - the number of columns and their global indices
634 
635    Notes:
636    The user must allocate space (m*n Scalars) for the values, v.
637    The values, v, are then returned in a row-oriented format,
638    analogous to that used by default in MatSetValues().
639 
640    MatGetValues() uses 0-based row and column numbers in
641    Fortran as well as in C.
642 
643    MatGetValues() requires that the matrix has been assembled
644    with MatAssemblyBegin()/MatAssemblyEnd().  Thus, calls to
645    MatSetValues() and MatGetValues() CANNOT be made in succession
646    without intermediate matrix assembly.
647 
648    Level: advanced
649 
650    Concepts: matrices^accessing values
651 
652 .seealso: MatGetRow(), MatGetSubmatrices(), MatSetValues()
653 @*/
654 int MatGetValues(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v)
655 {
656   int ierr;
657 
658   PetscFunctionBegin;
659   PetscValidHeaderSpecific(mat,MAT_COOKIE);
660   PetscValidType(mat);
661   MatPreallocated(mat);
662   PetscValidIntPointer(idxm);
663   PetscValidIntPointer(idxn);
664   PetscValidScalarPointer(v);
665   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
666   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
667   if (!mat->ops->getvalues) SETERRQ(PETSC_ERR_SUP,"");
668 
669   ierr = PLogEventBegin(MAT_GetValues,mat,0,0,0);CHKERRQ(ierr);
670   ierr = (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);CHKERRQ(ierr);
671   ierr = PLogEventEnd(MAT_GetValues,mat,0,0,0);CHKERRQ(ierr);
672   PetscFunctionReturn(0);
673 }
674 
675 #undef __FUNC__
676 #define __FUNC__ /*<a name=""></a>*/"MatSetLocalToGlobalMapping"
677 /*@
678    MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
679    the routine MatSetValuesLocal() to allow users to insert matrix entries
680    using a local (per-processor) numbering.
681 
682    Not Collective
683 
684    Input Parameters:
685 +  x - the matrix
686 -  mapping - mapping created with ISLocalToGlobalMappingCreate()
687              or ISLocalToGlobalMappingCreateIS()
688 
689    Level: intermediate
690 
691    Concepts: matrices^local to global mapping
692    Concepts: local to global mapping^for matrices
693 
694 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
695 @*/
696 int MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping)
697 {
698   int ierr;
699   PetscFunctionBegin;
700   PetscValidHeaderSpecific(x,MAT_COOKIE);
701   PetscValidType(x);
702   MatPreallocated(x);
703   PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE);
704   if (x->mapping) {
705     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
706   }
707 
708   if (x->ops->setlocaltoglobalmapping) {
709     ierr = (*x->ops->setlocaltoglobalmapping)(x,mapping);CHKERRQ(ierr);
710   } else {
711     x->mapping = mapping;
712     ierr = PetscObjectReference((PetscObject)mapping);CHKERRQ(ierr);
713   }
714   PetscFunctionReturn(0);
715 }
716 
717 #undef __FUNC__
718 #define __FUNC__ /*<a name=""></a>*/"MatSetLocalToGlobalMappingBlock"
719 /*@
720    MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use
721    by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
722    entries using a local (per-processor) numbering.
723 
724    Not Collective
725 
726    Input Parameters:
727 +  x - the matrix
728 -  mapping - mapping created with ISLocalToGlobalMappingCreate() or
729              ISLocalToGlobalMappingCreateIS()
730 
731    Level: intermediate
732 
733    Concepts: matrices^local to global mapping blocked
734    Concepts: local to global mapping^for matrices, blocked
735 
736 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
737            MatSetValuesBlocked(), MatSetValuesLocal()
738 @*/
739 int MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping mapping)
740 {
741   int ierr;
742   PetscFunctionBegin;
743   PetscValidHeaderSpecific(x,MAT_COOKIE);
744   PetscValidType(x);
745   MatPreallocated(x);
746   PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE);
747   if (x->bmapping) {
748     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
749   }
750 
751   x->bmapping = mapping;
752   ierr = PetscObjectReference((PetscObject)mapping);CHKERRQ(ierr);
753   PetscFunctionReturn(0);
754 }
755 
756 #undef __FUNC__
757 #define __FUNC__ /*<a name=""></a>*/"MatSetValuesLocal"
758 /*@
759    MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
760    using a local ordering of the nodes.
761 
762    Not Collective
763 
764    Input Parameters:
765 +  x - the matrix
766 .  nrow, irow - number of rows and their local indices
767 .  ncol, icol - number of columns and their local indices
768 .  y -  a logically two-dimensional array of values
769 -  addv - either INSERT_VALUES or ADD_VALUES, where
770    ADD_VALUES adds values to any existing entries, and
771    INSERT_VALUES replaces existing entries with new values
772 
773    Notes:
774    Before calling MatSetValuesLocal(), the user must first set the
775    local-to-global mapping by calling MatSetLocalToGlobalMapping().
776 
777    Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES
778    options cannot be mixed without intervening calls to the assembly
779    routines.
780 
781    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
782    MUST be called after all calls to MatSetValuesLocal() have been completed.
783 
784    Level: intermediate
785 
786    Concepts: matrices^putting entries in with local numbering
787 
788 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping()
789 @*/
790 int MatSetValuesLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,Scalar *y,InsertMode addv)
791 {
792   int ierr,irowm[2048],icolm[2048];
793 
794   PetscFunctionBegin;
795   PetscValidHeaderSpecific(mat,MAT_COOKIE);
796   PetscValidType(mat);
797   MatPreallocated(mat);
798   PetscValidIntPointer(irow);
799   PetscValidIntPointer(icol);
800   PetscValidScalarPointer(y);
801 
802   if (mat->insertmode == NOT_SET_VALUES) {
803     mat->insertmode = addv;
804   }
805 #if defined(PETSC_USE_BOPT_g)
806   else if (mat->insertmode != addv) {
807     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
808   }
809   if (!mat->ops->setvalueslocal && (nrow > 2048 || ncol > 2048)) {
810     SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
811   }
812   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
813 #endif
814 
815   if (mat->assembled) {
816     mat->was_assembled = PETSC_TRUE;
817     mat->assembled     = PETSC_FALSE;
818   }
819   ierr = PLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr);
820   if (!mat->ops->setvalueslocal) {
821     ierr = ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm);CHKERRQ(ierr);
822     ierr = ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);CHKERRQ(ierr);
823     ierr = (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr);
824   } else {
825     ierr = (*mat->ops->setvalueslocal)(mat,nrow,irow,ncol,icol,y,addv);CHKERRQ(ierr);
826   }
827   ierr = PLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr);
828   PetscFunctionReturn(0);
829 }
830 
831 #undef __FUNC__
832 #define __FUNC__ /*<a name=""></a>*/"MatSetValuesBlockedLocal"
833 /*@
834    MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
835    using a local ordering of the nodes a block at a time.
836 
837    Not Collective
838 
839    Input Parameters:
840 +  x - the matrix
841 .  nrow, irow - number of rows and their local indices
842 .  ncol, icol - number of columns and their local indices
843 .  y -  a logically two-dimensional array of values
844 -  addv - either INSERT_VALUES or ADD_VALUES, where
845    ADD_VALUES adds values to any existing entries, and
846    INSERT_VALUES replaces existing entries with new values
847 
848    Notes:
849    Before calling MatSetValuesBlockedLocal(), the user must first set the
850    local-to-global mapping by calling MatSetLocalToGlobalMappingBlock(),
851    where the mapping MUST be set for matrix blocks, not for matrix elements.
852 
853    Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
854    options cannot be mixed without intervening calls to the assembly
855    routines.
856 
857    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
858    MUST be called after all calls to MatSetValuesBlockedLocal() have been completed.
859 
860    Level: intermediate
861 
862    Concepts: matrices^putting blocked values in with local numbering
863 
864 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked()
865 @*/
866 int MatSetValuesBlockedLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,Scalar *y,InsertMode addv)
867 {
868   int ierr,irowm[2048],icolm[2048];
869 
870   PetscFunctionBegin;
871   PetscValidHeaderSpecific(mat,MAT_COOKIE);
872   PetscValidType(mat);
873   MatPreallocated(mat);
874   PetscValidIntPointer(irow);
875   PetscValidIntPointer(icol);
876   PetscValidScalarPointer(y);
877   if (mat->insertmode == NOT_SET_VALUES) {
878     mat->insertmode = addv;
879   }
880 #if defined(PETSC_USE_BOPT_g)
881   else if (mat->insertmode != addv) {
882     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
883   }
884   if (!mat->bmapping) {
885     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with MatSetLocalToGlobalMappingBlock()");
886   }
887   if (nrow > 2048 || ncol > 2048) {
888     SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
889   }
890   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
891 #endif
892 
893   if (mat->assembled) {
894     mat->was_assembled = PETSC_TRUE;
895     mat->assembled     = PETSC_FALSE;
896   }
897   ierr = PLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr);
898   ierr = ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm);CHKERRQ(ierr);
899   ierr = ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm);CHKERRQ(ierr);
900   ierr = (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr);
901   ierr = PLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr);
902   PetscFunctionReturn(0);
903 }
904 
905 /* --------------------------------------------------------*/
906 #undef __FUNC__
907 #define __FUNC__ /*<a name=""></a>*/"MatMult"
908 /*@
909    MatMult - Computes the matrix-vector product, y = Ax.
910 
911    Collective on Mat and Vec
912 
913    Input Parameters:
914 +  mat - the matrix
915 -  x   - the vector to be multilplied
916 
917    Output Parameters:
918 .  y - the result
919 
920    Notes:
921    The vectors x and y cannot be the same.  I.e., one cannot
922    call MatMult(A,y,y).
923 
924    Level: beginner
925 
926    Concepts: matrix-vector product
927 
928 .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
929 @*/
930 int MatMult(Mat mat,Vec x,Vec y)
931 {
932   int ierr;
933 
934   PetscFunctionBegin;
935   PetscValidHeaderSpecific(mat,MAT_COOKIE);
936   PetscValidType(mat);
937   MatPreallocated(mat);
938   PetscValidHeaderSpecific(x,VEC_COOKIE);
939   PetscValidHeaderSpecific(y,VEC_COOKIE);
940 
941   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
942   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
943   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
944   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
945   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
946   if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %d %d",mat->m,y->n);
947 
948   if (mat->nullsp) {
949     ierr = MatNullSpaceRemove(mat->nullsp,x,&x);CHKERRQ(ierr);
950   }
951 
952   ierr = PLogEventBegin(MAT_Mult,mat,x,y,0);CHKERRQ(ierr);
953   ierr = (*mat->ops->mult)(mat,x,y);CHKERRQ(ierr);
954   ierr = PLogEventEnd(MAT_Mult,mat,x,y,0);CHKERRQ(ierr);
955 
956   if (mat->nullsp) {
957     ierr = MatNullSpaceRemove(mat->nullsp,y,PETSC_NULL);CHKERRQ(ierr);
958   }
959   PetscFunctionReturn(0);
960 }
961 
962 #undef __FUNC__
963 #define __FUNC__ /*<a name=""></a>*/"MatMultTranspose"
964 /*@
965    MatMultTranspose - Computes matrix transpose times a vector.
966 
967    Collective on Mat and Vec
968 
969    Input Parameters:
970 +  mat - the matrix
971 -  x   - the vector to be multilplied
972 
973    Output Parameters:
974 .  y - the result
975 
976    Notes:
977    The vectors x and y cannot be the same.  I.e., one cannot
978    call MatMultTranspose(A,y,y).
979 
980    Level: beginner
981 
982    Concepts: matrix vector product^transpose
983 
984 .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd()
985 @*/
986 int MatMultTranspose(Mat mat,Vec x,Vec y)
987 {
988   int ierr;
989 
990   PetscFunctionBegin;
991   PetscValidHeaderSpecific(mat,MAT_COOKIE);
992   PetscValidType(mat);
993   MatPreallocated(mat);
994   PetscValidHeaderSpecific(x,VEC_COOKIE);
995   PetscValidHeaderSpecific(y,VEC_COOKIE);
996 
997   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
998   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
999   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1000   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
1001   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
1002 
1003   ierr = PLogEventBegin(MAT_MultTranspose,mat,x,y,0);CHKERRQ(ierr);
1004   ierr = (*mat->ops->multtranspose)(mat,x,y);CHKERRQ(ierr);
1005   ierr = PLogEventEnd(MAT_MultTranspose,mat,x,y,0);CHKERRQ(ierr);
1006   PetscFunctionReturn(0);
1007 }
1008 
1009 #undef __FUNC__
1010 #define __FUNC__ /*<a name=""></a>*/"MatMultAdd"
1011 /*@
1012     MatMultAdd -  Computes v3 = v2 + A * v1.
1013 
1014     Collective on Mat and Vec
1015 
1016     Input Parameters:
1017 +   mat - the matrix
1018 -   v1, v2 - the vectors
1019 
1020     Output Parameters:
1021 .   v3 - the result
1022 
1023     Notes:
1024     The vectors v1 and v3 cannot be the same.  I.e., one cannot
1025     call MatMultAdd(A,v1,v2,v1).
1026 
1027     Level: beginner
1028 
1029     Concepts: matrix vector product^addition
1030 
1031 .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd()
1032 @*/
1033 int MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1034 {
1035   int ierr;
1036 
1037   PetscFunctionBegin;
1038   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1039   PetscValidType(mat);
1040   MatPreallocated(mat);
1041   PetscValidHeaderSpecific(v1,VEC_COOKIE);
1042   PetscValidHeaderSpecific(v2,VEC_COOKIE);
1043   PetscValidHeaderSpecific(v3,VEC_COOKIE);
1044 
1045   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1046   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1047   if (mat->N != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->N,v1->N);
1048   if (mat->M != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->M,v2->N);
1049   if (mat->M != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->M,v3->N);
1050   if (mat->m != v3->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: local dim %d %d",mat->m,v3->n);
1051   if (mat->m != v2->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: local dim %d %d",mat->m,v2->n);
1052   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
1053 
1054   ierr = PLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);CHKERRQ(ierr);
1055   ierr = (*mat->ops->multadd)(mat,v1,v2,v3);CHKERRQ(ierr);
1056   ierr = PLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);CHKERRQ(ierr);
1057   PetscFunctionReturn(0);
1058 }
1059 
1060 #undef __FUNC__
1061 #define __FUNC__ /*<a name=""></a>*/"MatMultTransposeAdd"
1062 /*@
1063    MatMultTransposeAdd - Computes v3 = v2 + A' * v1.
1064 
1065    Collective on Mat and Vec
1066 
1067    Input Parameters:
1068 +  mat - the matrix
1069 -  v1, v2 - the vectors
1070 
1071    Output Parameters:
1072 .  v3 - the result
1073 
1074    Notes:
1075    The vectors v1 and v3 cannot be the same.  I.e., one cannot
1076    call MatMultTransposeAdd(A,v1,v2,v1).
1077 
1078    Level: beginner
1079 
1080    Concepts: matrix vector product^transpose and addition
1081 
1082 .seealso: MatMultTranspose(), MatMultAdd(), MatMult()
1083 @*/
1084 int MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1085 {
1086   int ierr;
1087 
1088   PetscFunctionBegin;
1089   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1090   PetscValidType(mat);
1091   MatPreallocated(mat);
1092   PetscValidHeaderSpecific(v1,VEC_COOKIE);
1093   PetscValidHeaderSpecific(v2,VEC_COOKIE);
1094   PetscValidHeaderSpecific(v3,VEC_COOKIE);
1095 
1096   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1097   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1098   if (!mat->ops->multtransposeadd) SETERRQ(PETSC_ERR_SUP,"");
1099   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
1100   if (mat->M != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->M,v1->N);
1101   if (mat->N != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->N,v2->N);
1102   if (mat->N != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->N,v3->N);
1103 
1104   ierr = PLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);CHKERRQ(ierr);
1105   ierr = (*mat->ops->multtransposeadd)(mat,v1,v2,v3);CHKERRQ(ierr);
1106   ierr = PLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);CHKERRQ(ierr);
1107   PetscFunctionReturn(0);
1108 }
1109 /* ------------------------------------------------------------*/
1110 #undef __FUNC__
1111 #define __FUNC__ /*<a name=""></a>*/"MatGetInfo"
1112 /*@C
1113    MatGetInfo - Returns information about matrix storage (number of
1114    nonzeros, memory, etc.).
1115 
1116    Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used
1117    as the flag
1118 
1119    Input Parameters:
1120 .  mat - the matrix
1121 
1122    Output Parameters:
1123 +  flag - flag indicating the type of parameters to be returned
1124    (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
1125    MAT_GLOBAL_SUM - sum over all processors)
1126 -  info - matrix information context
1127 
1128    Notes:
1129    The MatInfo context contains a variety of matrix data, including
1130    number of nonzeros allocated and used, number of mallocs during
1131    matrix assembly, etc.  Additional information for factored matrices
1132    is provided (such as the fill ratio, number of mallocs during
1133    factorization, etc.).  Much of this info is printed to STDOUT
1134    when using the runtime options
1135 $       -log_info -mat_view_info
1136 
1137    Example for C/C++ Users:
1138    See the file ${PETSC_DIR}/include/petscmat.h for a complete list of
1139    data within the MatInfo context.  For example,
1140 .vb
1141       MatInfo info;
1142       Mat     A;
1143       double  mal, nz_a, nz_u;
1144 
1145       MatGetInfo(A,MAT_LOCAL,&info);
1146       mal  = info.mallocs;
1147       nz_a = info.nz_allocated;
1148 .ve
1149 
1150    Example for Fortran Users:
1151    Fortran users should declare info as a double precision
1152    array of dimension MAT_INFO_SIZE, and then extract the parameters
1153    of interest.  See the file ${PETSC_DIR}/include/finclude/petscmat.h
1154    a complete list of parameter names.
1155 .vb
1156       double  precision info(MAT_INFO_SIZE)
1157       double  precision mal, nz_a
1158       Mat     A
1159       integer ierr
1160 
1161       call MatGetInfo(A,MAT_LOCAL,info,ierr)
1162       mal = info(MAT_INFO_MALLOCS)
1163       nz_a = info(MAT_INFO_NZ_ALLOCATED)
1164 .ve
1165 
1166     Level: intermediate
1167 
1168     Concepts: matrices^getting information on
1169 
1170 @*/
1171 int MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
1172 {
1173   int ierr;
1174 
1175   PetscFunctionBegin;
1176   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1177   PetscValidType(mat);
1178   MatPreallocated(mat);
1179   PetscValidPointer(info);
1180   if (!mat->ops->getinfo) SETERRQ(PETSC_ERR_SUP,"");
1181   ierr = (*mat->ops->getinfo)(mat,flag,info);CHKERRQ(ierr);
1182   PetscFunctionReturn(0);
1183 }
1184 
1185 /* ----------------------------------------------------------*/
1186 #undef __FUNC__
1187 #define __FUNC__ /*<a name=""></a>*/"MatILUDTFactor"
1188 /*@C
1189    MatILUDTFactor - Performs a drop tolerance ILU factorization.
1190 
1191    Collective on Mat
1192 
1193    Input Parameters:
1194 +  mat - the matrix
1195 .  info - information about the factorization to be done
1196 .  row - row permutation
1197 -  col - column permutation
1198 
1199    Output Parameters:
1200 .  fact - the factored matrix
1201 
1202    Level: developer
1203 
1204    Notes:
1205    Most users should employ the simplified SLES interface for linear solvers
1206    instead of working directly with matrix algebra routines such as this.
1207    See, e.g., SLESCreate().
1208 
1209    This is currently only supported for the SeqAIJ matrix format using code
1210    from Yousef Saad's SPARSEKIT2  package (translated to C with f2c) and/or
1211    Matlab. SPARSEKIT2 is copyrighted by Yousef Saad with the GNU copyright
1212    and thus can be distributed with PETSc.
1213 
1214     Concepts: matrices^ILUDT factorization
1215 
1216 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
1217 @*/
1218 int MatILUDTFactor(Mat mat,MatILUInfo *info,IS row,IS col,Mat *fact)
1219 {
1220   int ierr;
1221 
1222   PetscFunctionBegin;
1223   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1224   PetscValidType(mat);
1225   MatPreallocated(mat);
1226   PetscValidPointer(fact);
1227   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1228   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1229   if (!mat->ops->iludtfactor) SETERRQ(PETSC_ERR_SUP,"");
1230 
1231   ierr = PLogEventBegin(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr);
1232   ierr = (*mat->ops->iludtfactor)(mat,info,row,col,fact);CHKERRQ(ierr);
1233   ierr = PLogEventEnd(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr);
1234 
1235   PetscFunctionReturn(0);
1236 }
1237 
1238 #undef __FUNC__
1239 #define __FUNC__ /*<a name=""></a>*/"MatLUFactor"
1240 /*@
1241    MatLUFactor - Performs in-place LU factorization of matrix.
1242 
1243    Collective on Mat
1244 
1245    Input Parameters:
1246 +  mat - the matrix
1247 .  row - row permutation
1248 .  col - column permutation
1249 -  info - options for factorization, includes
1250 $          fill - expected fill as ratio of original fill.
1251 $          dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
1252 $                   Run with the option -log_info to determine an optimal value to use
1253 
1254    Notes:
1255    Most users should employ the simplified SLES interface for linear solvers
1256    instead of working directly with matrix algebra routines such as this.
1257    See, e.g., SLESCreate().
1258 
1259    This changes the state of the matrix to a factored matrix; it cannot be used
1260    for example with MatSetValues() unless one first calls MatSetUnfactored().
1261 
1262    Level: developer
1263 
1264    Concepts: matrices^LU factorization
1265 
1266 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(),
1267           MatGetOrdering(), MatSetUnfactored()
1268 
1269 @*/
1270 int MatLUFactor(Mat mat,IS row,IS col,MatLUInfo *info)
1271 {
1272   int ierr;
1273 
1274   PetscFunctionBegin;
1275   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1276   PetscValidType(mat);
1277   MatPreallocated(mat);
1278   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1279   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1280   if (!mat->ops->lufactor) SETERRQ(PETSC_ERR_SUP,"");
1281 
1282   ierr = PLogEventBegin(MAT_LUFactor,mat,row,col,0);CHKERRQ(ierr);
1283   ierr = (*mat->ops->lufactor)(mat,row,col,info);CHKERRQ(ierr);
1284   ierr = PLogEventEnd(MAT_LUFactor,mat,row,col,0);CHKERRQ(ierr);
1285   PetscFunctionReturn(0);
1286 }
1287 
1288 #undef __FUNC__
1289 #define __FUNC__ /*<a name=""></a>*/"MatILUFactor"
1290 /*@
1291    MatILUFactor - Performs in-place ILU factorization of matrix.
1292 
1293    Collective on Mat
1294 
1295    Input Parameters:
1296 +  mat - the matrix
1297 .  row - row permutation
1298 .  col - column permutation
1299 -  info - structure containing
1300 $      levels - number of levels of fill.
1301 $      expected fill - as ratio of original fill.
1302 $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
1303                 missing diagonal entries)
1304 
1305    Notes:
1306    Probably really in-place only when level of fill is zero, otherwise allocates
1307    new space to store factored matrix and deletes previous memory.
1308 
1309    Most users should employ the simplified SLES interface for linear solvers
1310    instead of working directly with matrix algebra routines such as this.
1311    See, e.g., SLESCreate().
1312 
1313    Level: developer
1314 
1315    Concepts: matrices^ILU factorization
1316 
1317 .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
1318 @*/
1319 int MatILUFactor(Mat mat,IS row,IS col,MatILUInfo *info)
1320 {
1321   int ierr;
1322 
1323   PetscFunctionBegin;
1324   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1325   PetscValidType(mat);
1326   MatPreallocated(mat);
1327   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
1328   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1329   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1330   if (!mat->ops->ilufactor) SETERRQ(PETSC_ERR_SUP,"");
1331 
1332   ierr = PLogEventBegin(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr);
1333   ierr = (*mat->ops->ilufactor)(mat,row,col,info);CHKERRQ(ierr);
1334   ierr = PLogEventEnd(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr);
1335   PetscFunctionReturn(0);
1336 }
1337 
1338 #undef __FUNC__
1339 #define __FUNC__ /*<a name=""></a>*/"MatLUFactorSymbolic"
1340 /*@
1341    MatLUFactorSymbolic - Performs symbolic LU factorization of matrix.
1342    Call this routine before calling MatLUFactorNumeric().
1343 
1344    Collective on Mat
1345 
1346    Input Parameters:
1347 +  mat - the matrix
1348 .  row, col - row and column permutations
1349 -  info - options for factorization, includes
1350 $          fill - expected fill as ratio of original fill.
1351 $          dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
1352 $                   Run with the option -log_info to determine an optimal value to use
1353 
1354    Output Parameter:
1355 .  fact - new matrix that has been symbolically factored
1356 
1357    Notes:
1358    See the users manual for additional information about
1359    choosing the fill factor for better efficiency.
1360 
1361    Most users should employ the simplified SLES interface for linear solvers
1362    instead of working directly with matrix algebra routines such as this.
1363    See, e.g., SLESCreate().
1364 
1365    Level: developer
1366 
1367    Concepts: matrices^LU symbolic factorization
1368 
1369 .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor()
1370 @*/
1371 int MatLUFactorSymbolic(Mat mat,IS row,IS col,MatLUInfo *info,Mat *fact)
1372 {
1373   int ierr;
1374 
1375   PetscFunctionBegin;
1376   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1377   PetscValidType(mat);
1378   MatPreallocated(mat);
1379   PetscValidPointer(fact);
1380   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1381   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1382   if (!mat->ops->lufactorsymbolic) SETERRQ(PETSC_ERR_SUP,"");
1383 
1384   ierr = PLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr);
1385   ierr = (*mat->ops->lufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr);
1386   ierr = PLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr);
1387   PetscFunctionReturn(0);
1388 }
1389 
1390 #undef __FUNC__
1391 #define __FUNC__ /*<a name=""></a>*/"MatLUFactorNumeric"
1392 /*@
1393    MatLUFactorNumeric - Performs numeric LU factorization of a matrix.
1394    Call this routine after first calling MatLUFactorSymbolic().
1395 
1396    Collective on Mat
1397 
1398    Input Parameters:
1399 +  mat - the matrix
1400 -  fact - the matrix generated for the factor, from MatLUFactorSymbolic()
1401 
1402    Notes:
1403    See MatLUFactor() for in-place factorization.  See
1404    MatCholeskyFactorNumeric() for the symmetric, positive definite case.
1405 
1406    Most users should employ the simplified SLES interface for linear solvers
1407    instead of working directly with matrix algebra routines such as this.
1408    See, e.g., SLESCreate().
1409 
1410    Level: developer
1411 
1412    Concepts: matrices^LU numeric factorization
1413 
1414 .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
1415 @*/
1416 int MatLUFactorNumeric(Mat mat,Mat *fact)
1417 {
1418   int        ierr;
1419   PetscTruth flg;
1420 
1421   PetscFunctionBegin;
1422   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1423   PetscValidType(mat);
1424   MatPreallocated(mat);
1425   PetscValidPointer(fact);
1426   PetscValidHeaderSpecific(*fact,MAT_COOKIE);
1427   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1428   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1429     SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dimensions are different %d should = %d %d should = %d",
1430             mat->M,(*fact)->M,mat->N,(*fact)->N);
1431   }
1432   if (!(*fact)->ops->lufactornumeric) SETERRQ(PETSC_ERR_SUP,"");
1433 
1434   ierr = PLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr);
1435   ierr = (*(*fact)->ops->lufactornumeric)(mat,fact);CHKERRQ(ierr);
1436   ierr = PLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr);
1437   ierr = OptionsHasName(PETSC_NULL,"-mat_view_draw",&flg);CHKERRQ(ierr);
1438   if (flg) {
1439     ierr = OptionsHasName(PETSC_NULL,"-mat_view_contour",&flg);CHKERRQ(ierr);
1440     if (flg) {
1441       ierr = ViewerPushFormat(VIEWER_DRAW_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr);
1442     }
1443     ierr = MatView(*fact,VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
1444     ierr = ViewerFlush(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
1445     if (flg) {
1446       ierr = ViewerPopFormat(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
1447     }
1448   }
1449   PetscFunctionReturn(0);
1450 }
1451 
1452 #undef __FUNC__
1453 #define __FUNC__ /*<a name=""></a>*/"MatCholeskyFactor"
1454 /*@
1455    MatCholeskyFactor - Performs in-place Cholesky factorization of a
1456    symmetric matrix.
1457 
1458    Collective on Mat
1459 
1460    Input Parameters:
1461 +  mat - the matrix
1462 .  perm - row and column permutations
1463 -  f - expected fill as ratio of original fill
1464 
1465    Notes:
1466    See MatLUFactor() for the nonsymmetric case.  See also
1467    MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().
1468 
1469    Most users should employ the simplified SLES interface for linear solvers
1470    instead of working directly with matrix algebra routines such as this.
1471    See, e.g., SLESCreate().
1472 
1473    Level: developer
1474 
1475    Concepts: matrices^Cholesky factorization
1476 
1477 .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
1478           MatGetOrdering()
1479 
1480 @*/
1481 int MatCholeskyFactor(Mat mat,IS perm,PetscReal f)
1482 {
1483   int ierr;
1484 
1485   PetscFunctionBegin;
1486   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1487   PetscValidType(mat);
1488   MatPreallocated(mat);
1489   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
1490   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1491   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1492   if (!mat->ops->choleskyfactor) SETERRQ(PETSC_ERR_SUP,"");
1493 
1494   ierr = PLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);CHKERRQ(ierr);
1495   ierr = (*mat->ops->choleskyfactor)(mat,perm,f);CHKERRQ(ierr);
1496   ierr = PLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);CHKERRQ(ierr);
1497   PetscFunctionReturn(0);
1498 }
1499 
1500 #undef __FUNC__
1501 #define __FUNC__ /*<a name=""></a>*/"MatCholeskyFactorSymbolic"
1502 /*@
1503    MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
1504    of a symmetric matrix.
1505 
1506    Collective on Mat
1507 
1508    Input Parameters:
1509 +  mat - the matrix
1510 .  perm - row and column permutations
1511 -  f - expected fill as ratio of original
1512 
1513    Output Parameter:
1514 .  fact - the factored matrix
1515 
1516    Notes:
1517    See MatLUFactorSymbolic() for the nonsymmetric case.  See also
1518    MatCholeskyFactor() and MatCholeskyFactorNumeric().
1519 
1520    Most users should employ the simplified SLES interface for linear solvers
1521    instead of working directly with matrix algebra routines such as this.
1522    See, e.g., SLESCreate().
1523 
1524    Level: developer
1525 
1526    Concepts: matrices^Cholesky symbolic factorization
1527 
1528 .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
1529           MatGetOrdering()
1530 
1531 @*/
1532 int MatCholeskyFactorSymbolic(Mat mat,IS perm,PetscReal f,Mat *fact)
1533 {
1534   int ierr;
1535 
1536   PetscFunctionBegin;
1537   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1538   PetscValidType(mat);
1539   MatPreallocated(mat);
1540   PetscValidPointer(fact);
1541   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
1542   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1543   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1544   if (!mat->ops->choleskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,"");
1545 
1546   ierr = PLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr);
1547   ierr = (*mat->ops->choleskyfactorsymbolic)(mat,perm,f,fact);CHKERRQ(ierr);
1548   ierr = PLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr);
1549   PetscFunctionReturn(0);
1550 }
1551 
1552 #undef __FUNC__
1553 #define __FUNC__ /*<a name=""></a>*/"MatCholeskyFactorNumeric"
1554 /*@
1555    MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
1556    of a symmetric matrix. Call this routine after first calling
1557    MatCholeskyFactorSymbolic().
1558 
1559    Collective on Mat
1560 
1561    Input Parameter:
1562 .  mat - the initial matrix
1563 
1564    Output Parameter:
1565 .  fact - the factored matrix
1566 
1567    Notes:
1568    Most users should employ the simplified SLES interface for linear solvers
1569    instead of working directly with matrix algebra routines such as this.
1570    See, e.g., SLESCreate().
1571 
1572    Level: developer
1573 
1574    Concepts: matrices^Cholesky numeric factorization
1575 
1576 .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
1577 @*/
1578 int MatCholeskyFactorNumeric(Mat mat,Mat *fact)
1579 {
1580   int ierr;
1581 
1582   PetscFunctionBegin;
1583   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1584   PetscValidType(mat);
1585   MatPreallocated(mat);
1586   PetscValidPointer(fact);
1587   if (!mat->ops->choleskyfactornumeric) SETERRQ(PETSC_ERR_SUP,"");
1588   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1589   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1590     SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dim %d should = %d %d should = %d",
1591             mat->M,(*fact)->M,mat->N,(*fact)->N);
1592   }
1593 
1594   ierr = PLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr);
1595   ierr = (*mat->ops->choleskyfactornumeric)(mat,fact);CHKERRQ(ierr);
1596   ierr = PLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr);
1597   PetscFunctionReturn(0);
1598 }
1599 
1600 /* ----------------------------------------------------------------*/
1601 #undef __FUNC__
1602 #define __FUNC__ /*<a name=""></a>*/"MatSolve"
1603 /*@
1604    MatSolve - Solves A x = b, given a factored matrix.
1605 
1606    Collective on Mat and Vec
1607 
1608    Input Parameters:
1609 +  mat - the factored matrix
1610 -  b - the right-hand-side vector
1611 
1612    Output Parameter:
1613 .  x - the result vector
1614 
1615    Notes:
1616    The vectors b and x cannot be the same.  I.e., one cannot
1617    call MatSolve(A,x,x).
1618 
1619    Notes:
1620    Most users should employ the simplified SLES interface for linear solvers
1621    instead of working directly with matrix algebra routines such as this.
1622    See, e.g., SLESCreate().
1623 
1624    Level: developer
1625 
1626    Concepts: matrices^triangular solves
1627 
1628 .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd()
1629 @*/
1630 int MatSolve(Mat mat,Vec b,Vec x)
1631 {
1632   int ierr;
1633 
1634   PetscFunctionBegin;
1635   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1636   PetscValidType(mat);
1637   MatPreallocated(mat);
1638   PetscValidHeaderSpecific(b,VEC_COOKIE);
1639   PetscValidHeaderSpecific(x,VEC_COOKIE);
1640   PetscCheckSameComm(mat,b);
1641   PetscCheckSameComm(mat,x);
1642   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
1643   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
1644   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1645   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1646   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1647   if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0);
1648 
1649   if (!mat->ops->solve) SETERRQ(PETSC_ERR_SUP,"");
1650   ierr = PLogEventBegin(MAT_Solve,mat,b,x,0);CHKERRQ(ierr);
1651   ierr = (*mat->ops->solve)(mat,b,x);CHKERRQ(ierr);
1652   ierr = PLogEventEnd(MAT_Solve,mat,b,x,0);CHKERRQ(ierr);
1653   PetscFunctionReturn(0);
1654 }
1655 
1656 #undef __FUNC__
1657 #define __FUNC__ /*<a name=""></a>*/"MatForwardSolve"
1658 /* @
1659    MatForwardSolve - Solves L x = b, given a factored matrix, A = LU.
1660 
1661    Collective on Mat and Vec
1662 
1663    Input Parameters:
1664 +  mat - the factored matrix
1665 -  b - the right-hand-side vector
1666 
1667    Output Parameter:
1668 .  x - the result vector
1669 
1670    Notes:
1671    MatSolve() should be used for most applications, as it performs
1672    a forward solve followed by a backward solve.
1673 
1674    The vectors b and x cannot be the same.  I.e., one cannot
1675    call MatForwardSolve(A,x,x).
1676 
1677    Most users should employ the simplified SLES interface for linear solvers
1678    instead of working directly with matrix algebra routines such as this.
1679    See, e.g., SLESCreate().
1680 
1681    Level: developer
1682 
1683    Concepts: matrices^forward solves
1684 
1685 .seealso: MatSolve(), MatBackwardSolve()
1686 @ */
1687 int MatForwardSolve(Mat mat,Vec b,Vec x)
1688 {
1689   int ierr;
1690 
1691   PetscFunctionBegin;
1692   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1693   PetscValidType(mat);
1694   MatPreallocated(mat);
1695   PetscValidHeaderSpecific(b,VEC_COOKIE);
1696   PetscValidHeaderSpecific(x,VEC_COOKIE);
1697   PetscCheckSameComm(mat,b);
1698   PetscCheckSameComm(mat,x);
1699   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
1700   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
1701   if (!mat->ops->forwardsolve) SETERRQ(PETSC_ERR_SUP,"");
1702   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1703   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1704   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1705 
1706   ierr = PLogEventBegin(MAT_ForwardSolve,mat,b,x,0);CHKERRQ(ierr);
1707   ierr = (*mat->ops->forwardsolve)(mat,b,x);CHKERRQ(ierr);
1708   ierr = PLogEventEnd(MAT_ForwardSolve,mat,b,x,0);CHKERRQ(ierr);
1709   PetscFunctionReturn(0);
1710 }
1711 
1712 #undef __FUNC__
1713 #define __FUNC__ /*<a name=""></a>*/"MatBackwardSolve"
1714 /* @
1715    MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.
1716 
1717    Collective on Mat and Vec
1718 
1719    Input Parameters:
1720 +  mat - the factored matrix
1721 -  b - the right-hand-side vector
1722 
1723    Output Parameter:
1724 .  x - the result vector
1725 
1726    Notes:
1727    MatSolve() should be used for most applications, as it performs
1728    a forward solve followed by a backward solve.
1729 
1730    The vectors b and x cannot be the same.  I.e., one cannot
1731    call MatBackwardSolve(A,x,x).
1732 
1733    Most users should employ the simplified SLES interface for linear solvers
1734    instead of working directly with matrix algebra routines such as this.
1735    See, e.g., SLESCreate().
1736 
1737    Level: developer
1738 
1739    Concepts: matrices^backward solves
1740 
1741 .seealso: MatSolve(), MatForwardSolve()
1742 @ */
1743 int MatBackwardSolve(Mat mat,Vec b,Vec x)
1744 {
1745   int ierr;
1746 
1747   PetscFunctionBegin;
1748   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1749   PetscValidType(mat);
1750   MatPreallocated(mat);
1751   PetscValidHeaderSpecific(b,VEC_COOKIE);
1752   PetscValidHeaderSpecific(x,VEC_COOKIE);
1753   PetscCheckSameComm(mat,b);
1754   PetscCheckSameComm(mat,x);
1755   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
1756   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
1757   if (!mat->ops->backwardsolve) SETERRQ(PETSC_ERR_SUP,"");
1758   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1759   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1760   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1761 
1762   ierr = PLogEventBegin(MAT_BackwardSolve,mat,b,x,0);CHKERRQ(ierr);
1763   ierr = (*mat->ops->backwardsolve)(mat,b,x);CHKERRQ(ierr);
1764   ierr = PLogEventEnd(MAT_BackwardSolve,mat,b,x,0);CHKERRQ(ierr);
1765   PetscFunctionReturn(0);
1766 }
1767 
1768 #undef __FUNC__
1769 #define __FUNC__ /*<a name=""></a>*/"MatSolveAdd"
1770 /*@
1771    MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.
1772 
1773    Collective on Mat and Vec
1774 
1775    Input Parameters:
1776 +  mat - the factored matrix
1777 .  b - the right-hand-side vector
1778 -  y - the vector to be added to
1779 
1780    Output Parameter:
1781 .  x - the result vector
1782 
1783    Notes:
1784    The vectors b and x cannot be the same.  I.e., one cannot
1785    call MatSolveAdd(A,x,y,x).
1786 
1787    Most users should employ the simplified SLES interface for linear solvers
1788    instead of working directly with matrix algebra routines such as this.
1789    See, e.g., SLESCreate().
1790 
1791    Level: developer
1792 
1793    Concepts: matrices^triangular solves
1794 
1795 .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd()
1796 @*/
1797 int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
1798 {
1799   Scalar one = 1.0;
1800   Vec    tmp;
1801   int    ierr;
1802 
1803   PetscFunctionBegin;
1804   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1805   PetscValidType(mat);
1806   MatPreallocated(mat);
1807   PetscValidHeaderSpecific(y,VEC_COOKIE);
1808   PetscValidHeaderSpecific(b,VEC_COOKIE);
1809   PetscValidHeaderSpecific(x,VEC_COOKIE);
1810   PetscCheckSameComm(mat,b);
1811   PetscCheckSameComm(mat,y);
1812   PetscCheckSameComm(mat,x);
1813   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
1814   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
1815   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1816   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1817   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
1818   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1819   if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n);
1820 
1821   ierr = PLogEventBegin(MAT_SolveAdd,mat,b,x,y);CHKERRQ(ierr);
1822   if (mat->ops->solveadd)  {
1823     ierr = (*mat->ops->solveadd)(mat,b,y,x);CHKERRQ(ierr);
1824   } else {
1825     /* do the solve then the add manually */
1826     if (x != y) {
1827       ierr = MatSolve(mat,b,x);CHKERRQ(ierr);
1828       ierr = VecAXPY(&one,y,x);CHKERRQ(ierr);
1829     } else {
1830       ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr);
1831       PLogObjectParent(mat,tmp);
1832       ierr = VecCopy(x,tmp);CHKERRQ(ierr);
1833       ierr = MatSolve(mat,b,x);CHKERRQ(ierr);
1834       ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr);
1835       ierr = VecDestroy(tmp);CHKERRQ(ierr);
1836     }
1837   }
1838   ierr = PLogEventEnd(MAT_SolveAdd,mat,b,x,y);CHKERRQ(ierr);
1839   PetscFunctionReturn(0);
1840 }
1841 
1842 #undef __FUNC__
1843 #define __FUNC__ /*<a name=""></a>*/"MatSolveTranspose"
1844 /*@
1845    MatSolveTranspose - Solves A' x = b, given a factored matrix.
1846 
1847    Collective on Mat and Vec
1848 
1849    Input Parameters:
1850 +  mat - the factored matrix
1851 -  b - the right-hand-side vector
1852 
1853    Output Parameter:
1854 .  x - the result vector
1855 
1856    Notes:
1857    The vectors b and x cannot be the same.  I.e., one cannot
1858    call MatSolveTranspose(A,x,x).
1859 
1860    Most users should employ the simplified SLES interface for linear solvers
1861    instead of working directly with matrix algebra routines such as this.
1862    See, e.g., SLESCreate().
1863 
1864    Level: developer
1865 
1866    Concepts: matrices^triangular solves
1867 
1868 .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd()
1869 @*/
1870 int MatSolveTranspose(Mat mat,Vec b,Vec x)
1871 {
1872   int ierr;
1873 
1874   PetscFunctionBegin;
1875   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1876   PetscValidType(mat);
1877   MatPreallocated(mat);
1878   PetscValidHeaderSpecific(b,VEC_COOKIE);
1879   PetscValidHeaderSpecific(x,VEC_COOKIE);
1880   PetscCheckSameComm(mat,b);
1881   PetscCheckSameComm(mat,x);
1882   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
1883   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
1884   if (!mat->ops->solvetranspose) SETERRQ(PETSC_ERR_SUP,"");
1885   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
1886   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);
1887 
1888   ierr = PLogEventBegin(MAT_SolveTranspose,mat,b,x,0);CHKERRQ(ierr);
1889   ierr = (*mat->ops->solvetranspose)(mat,b,x);CHKERRQ(ierr);
1890   ierr = PLogEventEnd(MAT_SolveTranspose,mat,b,x,0);CHKERRQ(ierr);
1891   PetscFunctionReturn(0);
1892 }
1893 
1894 #undef __FUNC__
1895 #define __FUNC__ /*<a name=""></a>*/"MatSolveTransposeAdd"
1896 /*@
1897    MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a
1898                       factored matrix.
1899 
1900    Collective on Mat and Vec
1901 
1902    Input Parameters:
1903 +  mat - the factored matrix
1904 .  b - the right-hand-side vector
1905 -  y - the vector to be added to
1906 
1907    Output Parameter:
1908 .  x - the result vector
1909 
1910    Notes:
1911    The vectors b and x cannot be the same.  I.e., one cannot
1912    call MatSolveTransposeAdd(A,x,y,x).
1913 
1914    Most users should employ the simplified SLES interface for linear solvers
1915    instead of working directly with matrix algebra routines such as this.
1916    See, e.g., SLESCreate().
1917 
1918    Level: developer
1919 
1920    Concepts: matrices^triangular solves
1921 
1922 .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose()
1923 @*/
1924 int MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x)
1925 {
1926   Scalar one = 1.0;
1927   int    ierr;
1928   Vec    tmp;
1929 
1930   PetscFunctionBegin;
1931   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1932   PetscValidType(mat);
1933   MatPreallocated(mat);
1934   PetscValidHeaderSpecific(y,VEC_COOKIE);
1935   PetscValidHeaderSpecific(b,VEC_COOKIE);
1936   PetscValidHeaderSpecific(x,VEC_COOKIE);
1937   PetscCheckSameComm(mat,b);
1938   PetscCheckSameComm(mat,y);
1939   PetscCheckSameComm(mat,x);
1940   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
1941   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
1942   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
1943   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);
1944   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
1945   if (x->n != y->n)   SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n);
1946 
1947   ierr = PLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);CHKERRQ(ierr);
1948   if (mat->ops->solvetransposeadd) {
1949     ierr = (*mat->ops->solvetransposeadd)(mat,b,y,x);CHKERRQ(ierr);
1950   } else {
1951     /* do the solve then the add manually */
1952     if (x != y) {
1953       ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr);
1954       ierr = VecAXPY(&one,y,x);CHKERRQ(ierr);
1955     } else {
1956       ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr);
1957       PLogObjectParent(mat,tmp);
1958       ierr = VecCopy(x,tmp);CHKERRQ(ierr);
1959       ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr);
1960       ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr);
1961       ierr = VecDestroy(tmp);CHKERRQ(ierr);
1962     }
1963   }
1964   ierr = PLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);CHKERRQ(ierr);
1965   PetscFunctionReturn(0);
1966 }
1967 /* ----------------------------------------------------------------*/
1968 
1969 #undef __FUNC__
1970 #define __FUNC__ /*<a name=""></a>*/"MatRelax"
1971 /*@
1972    MatRelax - Computes one relaxation sweep.
1973 
1974    Collective on Mat and Vec
1975 
1976    Input Parameters:
1977 +  mat - the matrix
1978 .  b - the right hand side
1979 .  omega - the relaxation factor
1980 .  flag - flag indicating the type of SOR (see below)
1981 .  shift -  diagonal shift
1982 -  its - the number of iterations
1983 
1984    Output Parameters:
1985 .  x - the solution (can contain an initial guess)
1986 
1987    SOR Flags:
1988 .     SOR_FORWARD_SWEEP - forward SOR
1989 .     SOR_BACKWARD_SWEEP - backward SOR
1990 .     SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
1991 .     SOR_LOCAL_FORWARD_SWEEP - local forward SOR
1992 .     SOR_LOCAL_BACKWARD_SWEEP - local forward SOR
1993 .     SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
1994 .     SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies
1995          upper/lower triangular part of matrix to
1996          vector (with omega)
1997 .     SOR_ZERO_INITIAL_GUESS - zero initial guess
1998 
1999    Notes:
2000    SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
2001    SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings
2002    on each processor.
2003 
2004    Application programmers will not generally use MatRelax() directly,
2005    but instead will employ the SLES/PC interface.
2006 
2007    Notes for Advanced Users:
2008    The flags are implemented as bitwise inclusive or operations.
2009    For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
2010    to specify a zero initial guess for SSOR.
2011 
2012    Most users should employ the simplified SLES interface for linear solvers
2013    instead of working directly with matrix algebra routines such as this.
2014    See, e.g., SLESCreate().
2015 
2016    Level: developer
2017 
2018    Concepts: matrices^relaxation
2019    Concepts: matrices^SOR
2020    Concepts: matrices^Gauss-Seidel
2021 
2022 @*/
2023 int MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,int its,Vec x)
2024 {
2025   int ierr;
2026 
2027   PetscFunctionBegin;
2028   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2029   PetscValidType(mat);
2030   MatPreallocated(mat);
2031   PetscValidHeaderSpecific(b,VEC_COOKIE);
2032   PetscValidHeaderSpecific(x,VEC_COOKIE);
2033   PetscCheckSameComm(mat,b);
2034   PetscCheckSameComm(mat,x);
2035   if (!mat->ops->relax) SETERRQ(PETSC_ERR_SUP,"");
2036   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2037   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2038   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
2039   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
2040   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
2041 
2042   ierr = PLogEventBegin(MAT_Relax,mat,b,x,0);CHKERRQ(ierr);
2043   ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,x);CHKERRQ(ierr);
2044   ierr = PLogEventEnd(MAT_Relax,mat,b,x,0);CHKERRQ(ierr);
2045   PetscFunctionReturn(0);
2046 }
2047 
2048 #undef __FUNC__
2049 #define __FUNC__ /*<a name=""></a>*/"MatCopy_Basic"
2050 /*
2051       Default matrix copy routine.
2052 */
2053 int MatCopy_Basic(Mat A,Mat B,MatStructure str)
2054 {
2055   int    ierr,i,rstart,rend,nz,*cwork;
2056   Scalar *vwork;
2057 
2058   PetscFunctionBegin;
2059   ierr = MatZeroEntries(B);CHKERRQ(ierr);
2060   ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
2061   for (i=rstart; i<rend; i++) {
2062     ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr);
2063     ierr = MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);CHKERRQ(ierr);
2064     ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr);
2065   }
2066   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2067   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2068   PetscFunctionReturn(0);
2069 }
2070 
2071 #undef __FUNC__
2072 #define __FUNC__ /*<a name=""></a>*/"MatCopy"
2073 /*@C
2074    MatCopy - Copys a matrix to another matrix.
2075 
2076    Collective on Mat
2077 
2078    Input Parameters:
2079 +  A - the matrix
2080 -  str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN
2081 
2082    Output Parameter:
2083 .  B - where the copy is put
2084 
2085    Notes:
2086    If you use SAME_NONZERO_PATTERN then the zero matrices had better have the
2087    same nonzero pattern or the routine will crash.
2088 
2089    MatCopy() copies the matrix entries of a matrix to another existing
2090    matrix (after first zeroing the second matrix).  A related routine is
2091    MatConvert(), which first creates a new matrix and then copies the data.
2092 
2093    Level: intermediate
2094 
2095    Concepts: matrices^copying
2096 
2097 .seealso: MatConvert()
2098 @*/
2099 int MatCopy(Mat A,Mat B,MatStructure str)
2100 {
2101   int ierr;
2102 
2103   PetscFunctionBegin;
2104   PetscValidHeaderSpecific(A,MAT_COOKIE);
2105   PetscValidHeaderSpecific(B,MAT_COOKIE);
2106   PetscValidType(A);
2107   MatPreallocated(A);
2108   PetscValidType(B);
2109   MatPreallocated(B);
2110   PetscCheckSameComm(A,B);
2111   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2112   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2113   if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim %d %d",A->M,B->M,
2114                                              A->N,B->N);
2115 
2116   ierr = PLogEventBegin(MAT_Copy,A,B,0,0);CHKERRQ(ierr);
2117   if (A->ops->copy) {
2118     ierr = (*A->ops->copy)(A,B,str);CHKERRQ(ierr);
2119   } else { /* generic conversion */
2120     ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr);
2121   }
2122   ierr = PLogEventEnd(MAT_Copy,A,B,0,0);CHKERRQ(ierr);
2123   PetscFunctionReturn(0);
2124 }
2125 
2126 #include "petscsys.h"
2127 PetscTruth MatConvertRegisterAllCalled = PETSC_FALSE;
2128 FList      MatConvertList              = 0;
2129 
2130 #undef __FUNC__
2131 #define __FUNC__ /*<a name=""></a>*/"MatConvertRegister"
2132 /*@C
2133     MatConvertRegister - Allows one to register a routine that reads matrices
2134         from a binary file for a particular matrix type.
2135 
2136   Not Collective
2137 
2138   Input Parameters:
2139 +   type - the type of matrix (defined in include/petscmat.h), for example, MATSEQAIJ.
2140 -   Converter - the function that reads the matrix from the binary file.
2141 
2142   Level: developer
2143 
2144 .seealso: MatConvertRegisterAll(), MatConvert()
2145 
2146 @*/
2147 int MatConvertRegister(char *sname,char *path,char *name,int (*function)(Mat,MatType,Mat*))
2148 {
2149   int  ierr;
2150   char fullname[256];
2151 
2152   PetscFunctionBegin;
2153   ierr = FListConcat(path,name,fullname);CHKERRQ(ierr);
2154   ierr = FListAdd(&MatConvertList,sname,fullname,(int (*)(void*))function);CHKERRQ(ierr);
2155   PetscFunctionReturn(0);
2156 }
2157 
2158 #undef __FUNC__
2159 #define __FUNC__ /*<a name=""></a>*/"MatConvert"
2160 /*@C
2161    MatConvert - Converts a matrix to another matrix, either of the same
2162    or different type.
2163 
2164    Collective on Mat
2165 
2166    Input Parameters:
2167 +  mat - the matrix
2168 -  newtype - new matrix type.  Use MATSAME to create a new matrix of the
2169    same type as the original matrix.
2170 
2171    Output Parameter:
2172 .  M - pointer to place new matrix
2173 
2174    Notes:
2175    MatConvert() first creates a new matrix and then copies the data from
2176    the first matrix.  A related routine is MatCopy(), which copies the matrix
2177    entries of one matrix to another already existing matrix context.
2178 
2179    Level: intermediate
2180 
2181    Concepts: matrices^converting between storage formats
2182 
2183 .seealso: MatCopy(), MatDuplicate()
2184 @*/
2185 int MatConvert(Mat mat,MatType newtype,Mat *M)
2186 {
2187   int        ierr;
2188   PetscTruth sametype,issame,flg;
2189   char       convname[256],mtype[256];
2190 
2191   PetscFunctionBegin;
2192   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2193   PetscValidType(mat);
2194   MatPreallocated(mat);
2195   PetscValidPointer(M);
2196   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2197   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2198 
2199   ierr = OptionsGetString(PETSC_NULL,"-matconvert_type",mtype,256,&flg);CHKERRQ(ierr);
2200   if (flg) {
2201     newtype = mtype;
2202   }
2203   ierr = PLogEventBegin(MAT_Convert,mat,0,0,0);CHKERRQ(ierr);
2204 
2205   ierr = PetscTypeCompare((PetscObject)mat,newtype,&sametype);CHKERRQ(ierr);
2206   ierr = PetscStrcmp(newtype,"same",&issame);CHKERRQ(ierr);
2207   if ((sametype || issame) && mat->ops->duplicate) {
2208     ierr = (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);CHKERRQ(ierr);
2209   } else {
2210     int (*conv)(Mat,MatType,Mat*);
2211     ierr = PetscStrcpy(convname,"MatConvertTo_");CHKERRQ(ierr);
2212     ierr = PetscStrcat(convname,newtype);CHKERRQ(ierr);
2213     ierr =  FListFind(mat->comm,MatConvertList,convname,(int(**)(void*))&conv);CHKERRQ(ierr);
2214     if (conv) {
2215       ierr = (*conv)(mat,newtype,M);CHKERRQ(ierr);
2216     } else {
2217       ierr = PetscStrcpy(convname,"MatConvert_");CHKERRQ(ierr);
2218       ierr = PetscStrcat(convname,mat->type_name);CHKERRQ(ierr);
2219       ierr = PetscStrcat(convname,"_");CHKERRQ(ierr);
2220       ierr = PetscStrcat(convname,newtype);CHKERRQ(ierr);
2221       ierr = PetscStrcat(convname,"_C");CHKERRQ(ierr);
2222       ierr = PetscObjectQueryFunction((PetscObject)mat,convname,(void**)&conv);CHKERRQ(ierr);
2223       if (conv) {
2224         ierr = (*conv)(mat,newtype,M);CHKERRQ(ierr);
2225       } else {
2226         if (mat->ops->convert) {
2227           ierr = (*mat->ops->convert)(mat,newtype,M);CHKERRQ(ierr);
2228         } else {
2229           ierr = MatConvert_Basic(mat,newtype,M);CHKERRQ(ierr);
2230         }
2231       }
2232     }
2233   }
2234   ierr = PLogEventEnd(MAT_Convert,mat,0,0,0);CHKERRQ(ierr);
2235   PetscFunctionReturn(0);
2236 }
2237 
2238 
2239 #undef __FUNC__
2240 #define __FUNC__ /*<a name=""></a>*/"MatDuplicate"
2241 /*@C
2242    MatDuplicate - Duplicates a matrix including the non-zero structure.
2243 
2244    Collective on Mat
2245 
2246    Input Parameters:
2247 +  mat - the matrix
2248 -  op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero
2249         values as well or not
2250 
2251    Output Parameter:
2252 .  M - pointer to place new matrix
2253 
2254    Level: intermediate
2255 
2256    Concepts: matrices^duplicating
2257 
2258 .seealso: MatCopy(), MatConvert()
2259 @*/
2260 int MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
2261 {
2262   int ierr;
2263 
2264   PetscFunctionBegin;
2265   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2266   PetscValidType(mat);
2267   MatPreallocated(mat);
2268   PetscValidPointer(M);
2269   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2270   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2271 
2272   *M  = 0;
2273   ierr = PLogEventBegin(MAT_Convert,mat,0,0,0);CHKERRQ(ierr);
2274   if (!mat->ops->duplicate) {
2275     SETERRQ(PETSC_ERR_SUP,"Not written for this matrix type");
2276   }
2277   ierr = (*mat->ops->duplicate)(mat,op,M);CHKERRQ(ierr);
2278   ierr = PLogEventEnd(MAT_Convert,mat,0,0,0);CHKERRQ(ierr);
2279   PetscFunctionReturn(0);
2280 }
2281 
2282 #undef __FUNC__
2283 #define __FUNC__ /*<a name=""></a>*/"MatGetDiagonal"
2284 /*@
2285    MatGetDiagonal - Gets the diagonal of a matrix.
2286 
2287    Collective on Mat and Vec
2288 
2289    Input Parameters:
2290 +  mat - the matrix
2291 -  v - the vector for storing the diagonal
2292 
2293    Output Parameter:
2294 .  v - the diagonal of the matrix
2295 
2296    Notes:
2297    For the SeqAIJ matrix format, this routine may also be called
2298    on a LU factored matrix; in that case it routines the reciprocal of
2299    the diagonal entries in U. It returns the entries permuted by the
2300    row and column permutation used during the symbolic factorization.
2301 
2302    Level: intermediate
2303 
2304    Concepts: matrices^accessing diagonals
2305 
2306 .seealso: MatGetRow(), MatGetSubmatrices(), MatGetSubmatrix(), MatGetRowMax()
2307 @*/
2308 int MatGetDiagonal(Mat mat,Vec v)
2309 {
2310   int ierr;
2311 
2312   PetscFunctionBegin;
2313   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2314   PetscValidType(mat);
2315   MatPreallocated(mat);
2316   PetscValidHeaderSpecific(v,VEC_COOKIE);
2317   /* PetscCheckSameComm(mat,v); Could be MPI vector but Seq matrix cause of two submatrix storage */
2318   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2319   if (!mat->ops->getdiagonal) SETERRQ(PETSC_ERR_SUP,"");
2320 
2321   ierr = (*mat->ops->getdiagonal)(mat,v);CHKERRQ(ierr);
2322   PetscFunctionReturn(0);
2323 }
2324 
2325 #undef __FUNC__
2326 #define __FUNC__ /*<a name=""></a>*/"MatGetRowMax"
2327 /*@
2328    MatGetRowMax - Gets the maximum value (in absolute value) of each
2329         row of the matrix
2330 
2331    Collective on Mat and Vec
2332 
2333    Input Parameters:
2334 .  mat - the matrix
2335 
2336    Output Parameter:
2337 .  v - the vector for storing the maximums
2338 
2339    Level: intermediate
2340 
2341    Concepts: matrices^getting row maximums
2342 
2343 .seealso: MatGetDiagonal(), MatGetSubmatrices(), MatGetSubmatrix()
2344 @*/
2345 int MatGetRowMax(Mat mat,Vec v)
2346 {
2347   int ierr;
2348 
2349   PetscFunctionBegin;
2350   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2351   PetscValidType(mat);
2352   MatPreallocated(mat);
2353   PetscValidHeaderSpecific(v,VEC_COOKIE);
2354   /* PetscCheckSameComm(mat,v); Could be MPI vector but Seq matrix cause of two submatrix storage */
2355   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2356   if (!mat->ops->getrowmax) SETERRQ(PETSC_ERR_SUP,"");
2357 
2358   ierr = (*mat->ops->getrowmax)(mat,v);CHKERRQ(ierr);
2359   PetscFunctionReturn(0);
2360 }
2361 
2362 #undef __FUNC__
2363 #define __FUNC__ /*<a name=""></a>*/"MatTranspose"
2364 /*@C
2365    MatTranspose - Computes an in-place or out-of-place transpose of a matrix.
2366 
2367    Collective on Mat
2368 
2369    Input Parameter:
2370 .  mat - the matrix to transpose
2371 
2372    Output Parameters:
2373 .  B - the transpose (or pass in PETSC_NULL for an in-place transpose)
2374 
2375    Level: intermediate
2376 
2377    Concepts: matrices^transposing
2378 
2379 .seealso: MatMultTranspose(), MatMultTransposeAdd()
2380 @*/
2381 int MatTranspose(Mat mat,Mat *B)
2382 {
2383   int ierr;
2384 
2385   PetscFunctionBegin;
2386   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2387   PetscValidType(mat);
2388   MatPreallocated(mat);
2389   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2390   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2391   if (!mat->ops->transpose) SETERRQ(PETSC_ERR_SUP,"");
2392   ierr = (*mat->ops->transpose)(mat,B);CHKERRQ(ierr);
2393   PetscFunctionReturn(0);
2394 }
2395 
2396 #undef __FUNC__
2397 #define __FUNC__ /*<a name=""></a>*/"MatPermute"
2398 /*@C
2399    MatPermute - Creates a new matrix with rows and columns permuted from the
2400    original.
2401 
2402    Collective on Mat
2403 
2404    Input Parameters:
2405 +  mat - the matrix to permute
2406 .  row - row permutation, each processor supplies only the permutation for its rows
2407 -  col - column permutation, each processor needs the entire column permutation, that is
2408          this is the same size as the total number of columns in the matrix
2409 
2410    Output Parameters:
2411 .  B - the permuted matrix
2412 
2413    Level: advanced
2414 
2415    Concepts: matrices^permuting
2416 
2417 .seealso: MatGetOrdering()
2418 @*/
2419 int MatPermute(Mat mat,IS row,IS col,Mat *B)
2420 {
2421   int ierr;
2422 
2423   PetscFunctionBegin;
2424   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2425   PetscValidType(mat);
2426   MatPreallocated(mat);
2427   PetscValidHeaderSpecific(row,IS_COOKIE);
2428   PetscValidHeaderSpecific(col,IS_COOKIE);
2429   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2430   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2431   if (!mat->ops->permute) SETERRQ(PETSC_ERR_SUP,"");
2432   ierr = (*mat->ops->permute)(mat,row,col,B);CHKERRQ(ierr);
2433   PetscFunctionReturn(0);
2434 }
2435 
2436 #undef __FUNC__
2437 #define __FUNC__ /*<a name=""></a>*/"MatEqual"
2438 /*@
2439    MatEqual - Compares two matrices.
2440 
2441    Collective on Mat
2442 
2443    Input Parameters:
2444 +  A - the first matrix
2445 -  B - the second matrix
2446 
2447    Output Parameter:
2448 .  flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.
2449 
2450    Level: intermediate
2451 
2452    Concepts: matrices^equality between
2453 @*/
2454 int MatEqual(Mat A,Mat B,PetscTruth *flg)
2455 {
2456   int ierr;
2457 
2458   PetscFunctionBegin;
2459   PetscValidHeaderSpecific(A,MAT_COOKIE);
2460   PetscValidHeaderSpecific(B,MAT_COOKIE);
2461   PetscValidType(A);
2462   MatPreallocated(A);
2463   PetscValidType(B);
2464   MatPreallocated(B);
2465   PetscValidIntPointer(flg);
2466   PetscCheckSameComm(A,B);
2467   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2468   if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2469   if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim %d %d %d %d",
2470                                              A->M,B->M,A->N,B->N);
2471   if (!A->ops->equal) SETERRQ(PETSC_ERR_SUP,"");
2472   ierr = (*A->ops->equal)(A,B,flg);CHKERRQ(ierr);
2473   PetscFunctionReturn(0);
2474 }
2475 
2476 #undef __FUNC__
2477 #define __FUNC__ /*<a name=""></a>*/"MatDiagonalScale"
2478 /*@
2479    MatDiagonalScale - Scales a matrix on the left and right by diagonal
2480    matrices that are stored as vectors.  Either of the two scaling
2481    matrices can be PETSC_NULL.
2482 
2483    Collective on Mat
2484 
2485    Input Parameters:
2486 +  mat - the matrix to be scaled
2487 .  l - the left scaling vector (or PETSC_NULL)
2488 -  r - the right scaling vector (or PETSC_NULL)
2489 
2490    Notes:
2491    MatDiagonalScale() computes A = LAR, where
2492    L = a diagonal matrix, R = a diagonal matrix
2493 
2494    Level: intermediate
2495 
2496    Concepts: matrices^diagonal scaling
2497    Concepts: diagonal scaling of matrices
2498 
2499 .seealso: MatScale()
2500 @*/
2501 int MatDiagonalScale(Mat mat,Vec l,Vec r)
2502 {
2503   int ierr;
2504 
2505   PetscFunctionBegin;
2506   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2507   PetscValidType(mat);
2508   MatPreallocated(mat);
2509   if (!mat->ops->diagonalscale) SETERRQ(PETSC_ERR_SUP,"");
2510   if (l) {PetscValidHeaderSpecific(l,VEC_COOKIE);PetscCheckSameComm(mat,l);}
2511   if (r) {PetscValidHeaderSpecific(r,VEC_COOKIE);PetscCheckSameComm(mat,r);}
2512   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2513   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2514 
2515   ierr = PLogEventBegin(MAT_Scale,mat,0,0,0);CHKERRQ(ierr);
2516   ierr = (*mat->ops->diagonalscale)(mat,l,r);CHKERRQ(ierr);
2517   ierr = PLogEventEnd(MAT_Scale,mat,0,0,0);CHKERRQ(ierr);
2518   PetscFunctionReturn(0);
2519 }
2520 
2521 #undef __FUNC__
2522 #define __FUNC__ /*<a name=""></a>*/"MatScale"
2523 /*@
2524     MatScale - Scales all elements of a matrix by a given number.
2525 
2526     Collective on Mat
2527 
2528     Input Parameters:
2529 +   mat - the matrix to be scaled
2530 -   a  - the scaling value
2531 
2532     Output Parameter:
2533 .   mat - the scaled matrix
2534 
2535     Level: intermediate
2536 
2537     Concepts: matrices^scaling all entries
2538 
2539 .seealso: MatDiagonalScale()
2540 @*/
2541 int MatScale(Scalar *a,Mat mat)
2542 {
2543   int ierr;
2544 
2545   PetscFunctionBegin;
2546   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2547   PetscValidType(mat);
2548   MatPreallocated(mat);
2549   PetscValidScalarPointer(a);
2550   if (!mat->ops->scale) SETERRQ(PETSC_ERR_SUP,"");
2551   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2552   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2553 
2554   ierr = PLogEventBegin(MAT_Scale,mat,0,0,0);CHKERRQ(ierr);
2555   ierr = (*mat->ops->scale)(a,mat);CHKERRQ(ierr);
2556   ierr = PLogEventEnd(MAT_Scale,mat,0,0,0);CHKERRQ(ierr);
2557   PetscFunctionReturn(0);
2558 }
2559 
2560 #undef __FUNC__
2561 #define __FUNC__ /*<a name=""></a>*/"MatNorm"
2562 /*@
2563    MatNorm - Calculates various norms of a matrix.
2564 
2565    Collective on Mat
2566 
2567    Input Parameters:
2568 +  mat - the matrix
2569 -  type - the type of norm, NORM_1, NORM_2, NORM_FROBENIUS, NORM_INFINITY
2570 
2571    Output Parameters:
2572 .  norm - the resulting norm
2573 
2574    Level: intermediate
2575 
2576    Concepts: matrices^norm
2577    Concepts: norm^of matrix
2578 @*/
2579 int MatNorm(Mat mat,NormType type,PetscReal *norm)
2580 {
2581   int ierr;
2582 
2583   PetscFunctionBegin;
2584   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2585   PetscValidType(mat);
2586   MatPreallocated(mat);
2587   PetscValidScalarPointer(norm);
2588 
2589   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2590   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2591   if (!mat->ops->norm) SETERRQ(PETSC_ERR_SUP,"Not for this matrix type");
2592   ierr = (*mat->ops->norm)(mat,type,norm);CHKERRQ(ierr);
2593   PetscFunctionReturn(0);
2594 }
2595 
2596 /*
2597      This variable is used to prevent counting of MatAssemblyBegin() that
2598    are called from within a MatAssemblyEnd().
2599 */
2600 static int MatAssemblyEnd_InUse = 0;
2601 #undef __FUNC__
2602 #define __FUNC__ /*<a name=""></a>*/"MatAssemblyBegin"
2603 /*@
2604    MatAssemblyBegin - Begins assembling the matrix.  This routine should
2605    be called after completing all calls to MatSetValues().
2606 
2607    Collective on Mat
2608 
2609    Input Parameters:
2610 +  mat - the matrix
2611 -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
2612 
2613    Notes:
2614    MatSetValues() generally caches the values.  The matrix is ready to
2615    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
2616    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
2617    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
2618    using the matrix.
2619 
2620    Level: beginner
2621 
2622    Concepts: matrices^assembling
2623 
2624 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled()
2625 @*/
2626 int MatAssemblyBegin(Mat mat,MatAssemblyType type)
2627 {
2628   int ierr;
2629 
2630   PetscFunctionBegin;
2631   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2632   PetscValidType(mat);
2633   MatPreallocated(mat);
2634   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.\nDid you forget to call MatSetUnfactored()?");
2635   if (mat->assembled) {
2636     mat->was_assembled = PETSC_TRUE;
2637     mat->assembled     = PETSC_FALSE;
2638   }
2639   if (!MatAssemblyEnd_InUse) {
2640     ierr = PLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);CHKERRQ(ierr);
2641     if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);}
2642     ierr = PLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);CHKERRQ(ierr);
2643   } else {
2644     if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);}
2645   }
2646   PetscFunctionReturn(0);
2647 }
2648 
2649 #undef __FUNC__
2650 #define __FUNC__ /*<a name=""></a>*/"MatAssembed"
2651 /*@
2652    MatAssembled - Indicates if a matrix has been assembled and is ready for
2653      use; for example, in matrix-vector product.
2654 
2655    Collective on Mat
2656 
2657    Input Parameter:
2658 .  mat - the matrix
2659 
2660    Output Parameter:
2661 .  assembled - PETSC_TRUE or PETSC_FALSE
2662 
2663    Level: advanced
2664 
2665    Concepts: matrices^assembled?
2666 
2667 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin()
2668 @*/
2669 int MatAssembled(Mat mat,PetscTruth *assembled)
2670 {
2671   PetscFunctionBegin;
2672   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2673   PetscValidType(mat);
2674   MatPreallocated(mat);
2675   *assembled = mat->assembled;
2676   PetscFunctionReturn(0);
2677 }
2678 
2679 #undef __FUNC__
2680 #define __FUNC__ /*<a name=""></a>*/"MatView_Private"
2681 /*
2682     Processes command line options to determine if/how a matrix
2683   is to be viewed. Called by MatAssemblyEnd() and MatLoad().
2684 */
2685 int MatView_Private(Mat mat)
2686 {
2687   int        ierr;
2688   PetscTruth flg;
2689 
2690   PetscFunctionBegin;
2691   ierr = OptionsHasName(mat->prefix,"-mat_view_info",&flg);CHKERRQ(ierr);
2692   if (flg) {
2693     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO,0);CHKERRQ(ierr);
2694     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2695     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2696   }
2697   ierr = OptionsHasName(mat->prefix,"-mat_view_info_detailed",&flg);CHKERRQ(ierr);
2698   if (flg) {
2699     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO_LONG,0);CHKERRQ(ierr);
2700     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2701     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2702   }
2703   ierr = OptionsHasName(mat->prefix,"-mat_view",&flg);CHKERRQ(ierr);
2704   if (flg) {
2705     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2706   }
2707   ierr = OptionsHasName(mat->prefix,"-mat_view_matlab",&flg);CHKERRQ(ierr);
2708   if (flg) {
2709     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_MATLAB,"M");CHKERRQ(ierr);
2710     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2711     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2712   }
2713   ierr = OptionsHasName(mat->prefix,"-mat_view_draw",&flg);CHKERRQ(ierr);
2714   if (flg) {
2715     ierr = OptionsHasName(mat->prefix,"-mat_view_contour",&flg);CHKERRQ(ierr);
2716     if (flg) {
2717       ViewerPushFormat(VIEWER_DRAW_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr);
2718     }
2719     ierr = MatView(mat,VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
2720     ierr = ViewerFlush(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
2721     if (flg) {
2722       ViewerPopFormat(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
2723     }
2724   }
2725   ierr = OptionsHasName(mat->prefix,"-mat_view_socket",&flg);CHKERRQ(ierr);
2726   if (flg) {
2727     ierr = MatView(mat,VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr);
2728     ierr = ViewerFlush(VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr);
2729   }
2730   ierr = OptionsHasName(mat->prefix,"-mat_view_binary",&flg);CHKERRQ(ierr);
2731   if (flg) {
2732     ierr = MatView(mat,VIEWER_BINARY_(mat->comm));CHKERRQ(ierr);
2733     ierr = ViewerFlush(VIEWER_BINARY_(mat->comm));CHKERRQ(ierr);
2734   }
2735   PetscFunctionReturn(0);
2736 }
2737 
2738 #undef __FUNC__
2739 #define __FUNC__ /*<a name=""></a>*/"MatAssemblyEnd"
2740 /*@
2741    MatAssemblyEnd - Completes assembling the matrix.  This routine should
2742    be called after MatAssemblyBegin().
2743 
2744    Collective on Mat
2745 
2746    Input Parameters:
2747 +  mat - the matrix
2748 -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
2749 
2750    Options Database Keys:
2751 +  -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
2752 .  -mat_view_info_detailed - Prints more detailed info
2753 .  -mat_view - Prints matrix in ASCII format
2754 .  -mat_view_matlab - Prints matrix in Matlab format
2755 .  -mat_view_draw - Draws nonzero structure of matrix, using MatView() and DrawOpenX().
2756 .  -display <name> - Sets display name (default is host)
2757 -  -draw_pause <sec> - Sets number of seconds to pause after display
2758 
2759    Notes:
2760    MatSetValues() generally caches the values.  The matrix is ready to
2761    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
2762    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
2763    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
2764    using the matrix.
2765 
2766    Level: beginner
2767 
2768 .seealso: MatAssemblyBegin(), MatSetValues(), DrawOpenX(), MatView(), MatAssembled()
2769 @*/
2770 int MatAssemblyEnd(Mat mat,MatAssemblyType type)
2771 {
2772   int        ierr;
2773   static int inassm = 0;
2774 
2775   PetscFunctionBegin;
2776   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2777   PetscValidType(mat);
2778   MatPreallocated(mat);
2779 
2780   inassm++;
2781   MatAssemblyEnd_InUse++;
2782   if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
2783     ierr = PLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);CHKERRQ(ierr);
2784     if (mat->ops->assemblyend) {
2785       ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr);
2786     }
2787     ierr = PLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);CHKERRQ(ierr);
2788   } else {
2789     if (mat->ops->assemblyend) {
2790       ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr);
2791     }
2792   }
2793 
2794   /* Flush assembly is not a true assembly */
2795   if (type != MAT_FLUSH_ASSEMBLY) {
2796     mat->assembled  = PETSC_TRUE; mat->num_ass++;
2797   }
2798   mat->insertmode = NOT_SET_VALUES;
2799   MatAssemblyEnd_InUse--;
2800 
2801   if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) {
2802     ierr = MatView_Private(mat);CHKERRQ(ierr);
2803   }
2804   inassm--;
2805   PetscFunctionReturn(0);
2806 }
2807 
2808 
2809 #undef __FUNC__
2810 #define __FUNC__ /*<a name=""></a>*/"MatCompress"
2811 /*@
2812    MatCompress - Tries to store the matrix in as little space as
2813    possible.  May fail if memory is already fully used, since it
2814    tries to allocate new space.
2815 
2816    Collective on Mat
2817 
2818    Input Parameters:
2819 .  mat - the matrix
2820 
2821    Level: advanced
2822 
2823 @*/
2824 int MatCompress(Mat mat)
2825 {
2826   int ierr;
2827 
2828   PetscFunctionBegin;
2829   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2830   PetscValidType(mat);
2831   MatPreallocated(mat);
2832   if (mat->ops->compress) {ierr = (*mat->ops->compress)(mat);CHKERRQ(ierr);}
2833   PetscFunctionReturn(0);
2834 }
2835 
2836 #undef __FUNC__
2837 #define __FUNC__ /*<a name=""></a>*/"MatSetOption"
2838 /*@
2839    MatSetOption - Sets a parameter option for a matrix. Some options
2840    may be specific to certain storage formats.  Some options
2841    determine how values will be inserted (or added). Sorted,
2842    row-oriented input will generally assemble the fastest. The default
2843    is row-oriented, nonsorted input.
2844 
2845    Collective on Mat
2846 
2847    Input Parameters:
2848 +  mat - the matrix
2849 -  option - the option, one of those listed below (and possibly others),
2850              e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR
2851 
2852    Options Describing Matrix Structure:
2853 +    MAT_SYMMETRIC - symmetric in terms of both structure and value
2854 -    MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure
2855 
2856    Options For Use with MatSetValues():
2857    Insert a logically dense subblock, which can be
2858 +    MAT_ROW_ORIENTED - row-oriented
2859 .    MAT_COLUMN_ORIENTED - column-oriented
2860 .    MAT_ROWS_SORTED - sorted by row
2861 .    MAT_ROWS_UNSORTED - not sorted by row
2862 .    MAT_COLUMNS_SORTED - sorted by column
2863 -    MAT_COLUMNS_UNSORTED - not sorted by column
2864 
2865    Not these options reflect the data you pass in with MatSetValues(); it has
2866    nothing to do with how the data is stored internally in the matrix
2867    data structure.
2868 
2869    When (re)assembling a matrix, we can restrict the input for
2870    efficiency/debugging purposes.  These options include
2871 +    MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be
2872         allowed if they generate a new nonzero
2873 .    MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed
2874 .    MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if
2875          they generate a nonzero in a new diagonal (for block diagonal format only)
2876 .    MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only)
2877 .    MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries
2878 .    MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry
2879 -    MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly
2880 
2881    Notes:
2882    Some options are relevant only for particular matrix types and
2883    are thus ignored by others.  Other options are not supported by
2884    certain matrix types and will generate an error message if set.
2885 
2886    If using a Fortran 77 module to compute a matrix, one may need to
2887    use the column-oriented option (or convert to the row-oriented
2888    format).
2889 
2890    MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion
2891    that would generate a new entry in the nonzero structure is instead
2892    ignored.  Thus, if memory has not alredy been allocated for this particular
2893    data, then the insertion is ignored. For dense matrices, in which
2894    the entire array is allocated, no entries are ever ignored.
2895 
2896    MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion
2897    that would generate a new entry in the nonzero structure instead produces
2898    an error. (Currently supported for AIJ and BAIJ formats only.)
2899    This is a useful flag when using SAME_NONZERO_PATTERN in calling
2900    SLESSetOperators() to ensure that the nonzero pattern truely does
2901    remain unchanged.
2902 
2903    MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion
2904    that would generate a new entry that has not been preallocated will
2905    instead produce an error. (Currently supported for AIJ and BAIJ formats
2906    only.) This is a useful flag when debugging matrix memory preallocation.
2907 
2908    MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for
2909    other processors should be dropped, rather than stashed.
2910    This is useful if you know that the "owning" processor is also
2911    always generating the correct matrix entries, so that PETSc need
2912    not transfer duplicate entries generated on another processor.
2913 
2914    MAT_USE_HASH_TABLE indicates that a hash table be used to improve the
2915    searches during matrix assembly. When this flag is set, the hash table
2916    is created during the first Matrix Assembly. This hash table is
2917    used the next time through, during MatSetVaules()/MatSetVaulesBlocked()
2918    to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag
2919    should be used with MAT_USE_HASH_TABLE flag. This option is currently
2920    supported by MATMPIBAIJ format only.
2921 
2922    MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries
2923    are kept in the nonzero structure
2924 
2925    MAT_IGNORE_ZERO_ENTRIES - when using ADD_VALUES for AIJ matrices this will stop
2926    zero values from creating a zero location in the matrix
2927 
2928    MAT_USE_INODES - indicates using inode version of the code - works with AIJ and
2929    ROWBS matrix types
2930 
2931    MAT_DO_NOT_USE_INODES - indicates not using inode version of the code - works
2932    with AIJ and ROWBS matrix types
2933 
2934    Level: intermediate
2935 
2936    Concepts: matrices^setting options
2937 
2938 @*/
2939 int MatSetOption(Mat mat,MatOption op)
2940 {
2941   int ierr;
2942 
2943   PetscFunctionBegin;
2944   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2945   PetscValidType(mat);
2946   MatPreallocated(mat);
2947   if (mat->ops->setoption) {ierr = (*mat->ops->setoption)(mat,op);CHKERRQ(ierr);}
2948   PetscFunctionReturn(0);
2949 }
2950 
2951 #undef __FUNC__
2952 #define __FUNC__ /*<a name=""></a>*/"MatZeroEntries"
2953 /*@
2954    MatZeroEntries - Zeros all entries of a matrix.  For sparse matrices
2955    this routine retains the old nonzero structure.
2956 
2957    Collective on Mat
2958 
2959    Input Parameters:
2960 .  mat - the matrix
2961 
2962    Level: intermediate
2963 
2964    Concepts: matrices^zeroing
2965 
2966 .seealso: MatZeroRows()
2967 @*/
2968 int MatZeroEntries(Mat mat)
2969 {
2970   int ierr;
2971 
2972   PetscFunctionBegin;
2973   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2974   PetscValidType(mat);
2975   MatPreallocated(mat);
2976   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2977   if (!mat->ops->zeroentries) SETERRQ(PETSC_ERR_SUP,"");
2978 
2979   ierr = PLogEventBegin(MAT_ZeroEntries,mat,0,0,0);CHKERRQ(ierr);
2980   ierr = (*mat->ops->zeroentries)(mat);CHKERRQ(ierr);
2981   ierr = PLogEventEnd(MAT_ZeroEntries,mat,0,0,0);CHKERRQ(ierr);
2982   PetscFunctionReturn(0);
2983 }
2984 
2985 #undef __FUNC__
2986 #define __FUNC__ /*<a name=""></a>*/"MatZeroRows"
2987 /*@C
2988    MatZeroRows - Zeros all entries (except possibly the main diagonal)
2989    of a set of rows of a matrix.
2990 
2991    Collective on Mat
2992 
2993    Input Parameters:
2994 +  mat - the matrix
2995 .  is - index set of rows to remove
2996 -  diag - pointer to value put in all diagonals of eliminated rows.
2997           Note that diag is not a pointer to an array, but merely a
2998           pointer to a single value.
2999 
3000    Notes:
3001    For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
3002    but does not release memory.  For the dense and block diagonal
3003    formats this does not alter the nonzero structure.
3004 
3005    If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
3006    of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
3007    merely zeroed.
3008 
3009    The user can set a value in the diagonal entry (or for the AIJ and
3010    row formats can optionally remove the main diagonal entry from the
3011    nonzero structure as well, by passing a null pointer (PETSC_NULL
3012    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
3013 
3014    For the parallel case, all processes that share the matrix (i.e.,
3015    those in the communicator used for matrix creation) MUST call this
3016    routine, regardless of whether any rows being zeroed are owned by
3017    them.
3018 
3019 
3020    Level: intermediate
3021 
3022    Concepts: matrices^zeroing rows
3023 
3024 .seealso: MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
3025 @*/
3026 int MatZeroRows(Mat mat,IS is,Scalar *diag)
3027 {
3028   int ierr;
3029 
3030   PetscFunctionBegin;
3031   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3032   PetscValidType(mat);
3033   MatPreallocated(mat);
3034   PetscValidHeaderSpecific(is,IS_COOKIE);
3035   if (diag) PetscValidScalarPointer(diag);
3036   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3037   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3038   if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,"");
3039 
3040   ierr = (*mat->ops->zerorows)(mat,is,diag);CHKERRQ(ierr);
3041   ierr = MatView_Private(mat);CHKERRQ(ierr);
3042   PetscFunctionReturn(0);
3043 }
3044 
3045 #undef __FUNC__
3046 #define __FUNC__ /*<a name=""></a>*/"MatZeroRowsLocal"
3047 /*@C
3048    MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
3049    of a set of rows of a matrix; using local numbering of rows.
3050 
3051    Collective on Mat
3052 
3053    Input Parameters:
3054 +  mat - the matrix
3055 .  is - index set of rows to remove
3056 -  diag - pointer to value put in all diagonals of eliminated rows.
3057           Note that diag is not a pointer to an array, but merely a
3058           pointer to a single value.
3059 
3060    Notes:
3061    Before calling MatZeroRowsLocal(), the user must first set the
3062    local-to-global mapping by calling MatSetLocalToGlobalMapping().
3063 
3064    For the AIJ matrix formats this removes the old nonzero structure,
3065    but does not release memory.  For the dense and block diagonal
3066    formats this does not alter the nonzero structure.
3067 
3068    The user can set a value in the diagonal entry (or for the AIJ and
3069    row formats can optionally remove the main diagonal entry from the
3070    nonzero structure as well, by passing a null pointer (PETSC_NULL
3071    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).
3072 
3073    Level: intermediate
3074 
3075    Concepts: matrices^zeroing
3076 
3077 .seealso: MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
3078 @*/
3079 int MatZeroRowsLocal(Mat mat,IS is,Scalar *diag)
3080 {
3081   int ierr;
3082   IS  newis;
3083 
3084   PetscFunctionBegin;
3085   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3086   PetscValidType(mat);
3087   MatPreallocated(mat);
3088   PetscValidHeaderSpecific(is,IS_COOKIE);
3089   if (diag) PetscValidScalarPointer(diag);
3090   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3091   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3092 
3093   if (mat->ops->zerorowslocal) {
3094     ierr = (*mat->ops->zerorowslocal)(mat,is,diag);CHKERRQ(ierr);
3095   } else {
3096     if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Need to provide local to global mapping to matrix first");
3097     ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);CHKERRQ(ierr);
3098     ierr = (*mat->ops->zerorows)(mat,newis,diag);CHKERRQ(ierr);
3099     ierr = ISDestroy(newis);CHKERRQ(ierr);
3100   }
3101   PetscFunctionReturn(0);
3102 }
3103 
3104 #undef __FUNC__
3105 #define __FUNC__ /*<a name=""></a>*/"MatGetSize"
3106 /*@
3107    MatGetSize - Returns the numbers of rows and columns in a matrix.
3108 
3109    Not Collective
3110 
3111    Input Parameter:
3112 .  mat - the matrix
3113 
3114    Output Parameters:
3115 +  m - the number of global rows
3116 -  n - the number of global columns
3117 
3118    Level: beginner
3119 
3120    Concepts: matrices^size
3121 
3122 .seealso: MatGetLocalSize()
3123 @*/
3124 int MatGetSize(Mat mat,int *m,int* n)
3125 {
3126   PetscFunctionBegin;
3127   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3128   if (m) *m = mat->M;
3129   if (n) *n = mat->N;
3130   PetscFunctionReturn(0);
3131 }
3132 
3133 #undef __FUNC__
3134 #define __FUNC__ /*<a name=""></a>*/"MatGetLocalSize"
3135 /*@
3136    MatGetLocalSize - Returns the number of rows and columns in a matrix
3137    stored locally.  This information may be implementation dependent, so
3138    use with care.
3139 
3140    Not Collective
3141 
3142    Input Parameters:
3143 .  mat - the matrix
3144 
3145    Output Parameters:
3146 +  m - the number of local rows
3147 -  n - the number of local columns
3148 
3149    Level: beginner
3150 
3151    Concepts: matrices^local size
3152 
3153 .seealso: MatGetSize()
3154 @*/
3155 int MatGetLocalSize(Mat mat,int *m,int* n)
3156 {
3157   PetscFunctionBegin;
3158   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3159   if (m) *m = mat->m;
3160   if (n) *n = mat->n;
3161   PetscFunctionReturn(0);
3162 }
3163 
3164 #undef __FUNC__
3165 #define __FUNC__ /*<a name=""></a>*/"MatGetOwnershipRange"
3166 /*@
3167    MatGetOwnershipRange - Returns the range of matrix rows owned by
3168    this processor, assuming that the matrix is laid out with the first
3169    n1 rows on the first processor, the next n2 rows on the second, etc.
3170    For certain parallel layouts this range may not be well defined.
3171 
3172    Not Collective
3173 
3174    Input Parameters:
3175 .  mat - the matrix
3176 
3177    Output Parameters:
3178 +  m - the global index of the first local row
3179 -  n - one more than the global index of the last local row
3180 
3181    Level: beginner
3182 
3183    Concepts: matrices^row ownership
3184 @*/
3185 int MatGetOwnershipRange(Mat mat,int *m,int* n)
3186 {
3187   int ierr;
3188 
3189   PetscFunctionBegin;
3190   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3191   PetscValidType(mat);
3192   MatPreallocated(mat);
3193   if (m) PetscValidIntPointer(m);
3194   if (n) PetscValidIntPointer(n);
3195   if (!mat->ops->getownershiprange) SETERRQ(PETSC_ERR_SUP,"");
3196   ierr = (*mat->ops->getownershiprange)(mat,m,n);CHKERRQ(ierr);
3197   PetscFunctionReturn(0);
3198 }
3199 
3200 #undef __FUNC__
3201 #define __FUNC__ /*<a name=""></a>*/"MatILUFactorSymbolic"
3202 /*@
3203    MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
3204    Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric()
3205    to complete the factorization.
3206 
3207    Collective on Mat
3208 
3209    Input Parameters:
3210 +  mat - the matrix
3211 .  row - row permutation
3212 .  column - column permutation
3213 -  info - structure containing
3214 $      levels - number of levels of fill.
3215 $      expected fill - as ratio of original fill.
3216 $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
3217                 missing diagonal entries)
3218 
3219    Output Parameters:
3220 .  fact - new matrix that has been symbolically factored
3221 
3222    Notes:
3223    See the users manual for additional information about
3224    choosing the fill factor for better efficiency.
3225 
3226    Most users should employ the simplified SLES interface for linear solvers
3227    instead of working directly with matrix algebra routines such as this.
3228    See, e.g., SLESCreate().
3229 
3230    Level: developer
3231 
3232   Concepts: matrices^symbolic LU factorization
3233   Concepts: matrices^factorization
3234   Concepts: LU^symbolic factorization
3235 
3236 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
3237           MatGetOrdering()
3238 
3239 @*/
3240 int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatILUInfo *info,Mat *fact)
3241 {
3242   int ierr;
3243 
3244   PetscFunctionBegin;
3245   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3246   PetscValidType(mat);
3247   MatPreallocated(mat);
3248   PetscValidPointer(fact);
3249   if (info && info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels of fill negative %d",info->levels);
3250   if (info && info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",info->fill);
3251   if (!mat->ops->ilufactorsymbolic) SETERRQ(PETSC_ERR_SUP,"Only MatCreateMPIRowbs() matrices support parallel ILU");
3252   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3253   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3254 
3255   ierr = PLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr);
3256   ierr = (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr);
3257   ierr = PLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr);
3258   PetscFunctionReturn(0);
3259 }
3260 
3261 #undef __FUNC__
3262 #define __FUNC__ /*<a name=""></a>*/"MatIncompleteCholeskyFactorSymbolic"
3263 /*@
3264    MatIncompleteCholeskyFactorSymbolic - Performs symbolic incomplete
3265    Cholesky factorization for a symmetric matrix.  Use
3266    MatCholeskyFactorNumeric() to complete the factorization.
3267 
3268    Collective on Mat
3269 
3270    Input Parameters:
3271 +  mat - the matrix
3272 .  perm - row and column permutation
3273 .  fill - levels of fill
3274 -  f - expected fill as ratio of original fill
3275 
3276    Output Parameter:
3277 .  fact - the factored matrix
3278 
3279    Notes:
3280    Currently only no-fill factorization is supported.
3281 
3282    Most users should employ the simplified SLES interface for linear solvers
3283    instead of working directly with matrix algebra routines such as this.
3284    See, e.g., SLESCreate().
3285 
3286    Level: developer
3287 
3288   Concepts: matrices^symbolic incomplete Cholesky factorization
3289   Concepts: matrices^factorization
3290   Concepts: Cholsky^symbolic factorization
3291 
3292 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor()
3293 @*/
3294 int MatIncompleteCholeskyFactorSymbolic(Mat mat,IS perm,PetscReal f,int fill,Mat *fact)
3295 {
3296   int ierr;
3297 
3298   PetscFunctionBegin;
3299   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3300   PetscValidType(mat);
3301   MatPreallocated(mat);
3302   PetscValidPointer(fact);
3303   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3304   if (fill < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Fill negative %d",fill);
3305   if (f < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",f);
3306   if (!mat->ops->incompletecholeskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,"Currently only MatCreateMPIRowbs() matrices support ICC in parallel");
3307   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3308 
3309   ierr = PLogEventBegin(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr);
3310   ierr = (*mat->ops->incompletecholeskyfactorsymbolic)(mat,perm,f,fill,fact);CHKERRQ(ierr);
3311   ierr = PLogEventEnd(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr);
3312   PetscFunctionReturn(0);
3313 }
3314 
3315 #undef __FUNC__
3316 #define __FUNC__ /*<a name=""></a>*/"MatGetArray"
3317 /*@C
3318    MatGetArray - Returns a pointer to the element values in the matrix.
3319    The result of this routine is dependent on the underlying matrix data
3320    structure, and may not even work for certain matrix types.  You MUST
3321    call MatRestoreArray() when you no longer need to access the array.
3322 
3323    Not Collective
3324 
3325    Input Parameter:
3326 .  mat - the matrix
3327 
3328    Output Parameter:
3329 .  v - the location of the values
3330 
3331    Currently returns an array only for the dense formats, giving access to
3332    the local portion of the matrix in the usual Fortran column-oriented format.
3333 
3334    Fortran Note:
3335    This routine is used differently from Fortran, e.g.,
3336 .vb
3337         Mat         mat
3338         Scalar      mat_array(1)
3339         PetscOffset i_mat
3340         int         ierr
3341         call MatGetArray(mat,mat_array,i_mat,ierr)
3342 
3343   C  Access first local entry in matrix; note that array is
3344   C  treated as one dimensional
3345         value = mat_array(i_mat + 1)
3346 
3347         [... other code ...]
3348         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3349 .ve
3350 
3351    See the Fortran chapter of the users manual and
3352    petsc/src/mat/examples/tests for details.
3353 
3354    Level: advanced
3355 
3356    Concepts: matrices^access array
3357 
3358 .seealso: MatRestoreArray(), MatGetArrayF90()
3359 @*/
3360 int MatGetArray(Mat mat,Scalar **v)
3361 {
3362   int ierr;
3363 
3364   PetscFunctionBegin;
3365   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3366   PetscValidType(mat);
3367   MatPreallocated(mat);
3368   PetscValidPointer(v);
3369   if (!mat->ops->getarray) SETERRQ(PETSC_ERR_SUP,"");
3370   ierr = (*mat->ops->getarray)(mat,v);CHKERRQ(ierr);
3371   PetscFunctionReturn(0);
3372 }
3373 
3374 #undef __FUNC__
3375 #define __FUNC__ /*<a name=""></a>*/"MatRestoreArray"
3376 /*@C
3377    MatRestoreArray - Restores the matrix after MatGetArray() has been called.
3378 
3379    Not Collective
3380 
3381    Input Parameter:
3382 +  mat - the matrix
3383 -  v - the location of the values
3384 
3385    Fortran Note:
3386    This routine is used differently from Fortran, e.g.,
3387 .vb
3388         Mat         mat
3389         Scalar      mat_array(1)
3390         PetscOffset i_mat
3391         int         ierr
3392         call MatGetArray(mat,mat_array,i_mat,ierr)
3393 
3394   C  Access first local entry in matrix; note that array is
3395   C  treated as one dimensional
3396         value = mat_array(i_mat + 1)
3397 
3398         [... other code ...]
3399         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3400 .ve
3401 
3402    See the Fortran chapter of the users manual and
3403    petsc/src/mat/examples/tests for details
3404 
3405    Level: advanced
3406 
3407 .seealso: MatGetArray(), MatRestoreArrayF90()
3408 @*/
3409 int MatRestoreArray(Mat mat,Scalar **v)
3410 {
3411   int ierr;
3412 
3413   PetscFunctionBegin;
3414   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3415   PetscValidType(mat);
3416   MatPreallocated(mat);
3417   PetscValidPointer(v);
3418   if (!mat->ops->restorearray) SETERRQ(PETSC_ERR_SUP,"");
3419   ierr = (*mat->ops->restorearray)(mat,v);CHKERRQ(ierr);
3420   PetscFunctionReturn(0);
3421 }
3422 
3423 #undef __FUNC__
3424 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrices"
3425 /*@C
3426    MatGetSubMatrices - Extracts several submatrices from a matrix. If submat
3427    points to an array of valid matrices, they may be reused to store the new
3428    submatrices.
3429 
3430    Collective on Mat
3431 
3432    Input Parameters:
3433 +  mat - the matrix
3434 .  n   - the number of submatrixes to be extracted (on this processor, may be zero)
3435 .  irow, icol - index sets of rows and columns to extract
3436 -  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3437 
3438    Output Parameter:
3439 .  submat - the array of submatrices
3440 
3441    Notes:
3442    MatGetSubMatrices() can extract only sequential submatrices
3443    (from both sequential and parallel matrices). Use MatGetSubMatrix()
3444    to extract a parallel submatrix.
3445 
3446    When extracting submatrices from a parallel matrix, each processor can
3447    form a different submatrix by setting the rows and columns of its
3448    individual index sets according to the local submatrix desired.
3449 
3450    When finished using the submatrices, the user should destroy
3451    them with MatDestroySubMatrices().
3452 
3453    MAT_REUSE_MATRIX can only be used when the nonzero structure of the
3454    original matrix has not changed from that last call to MatGetSubMatrices().
3455 
3456    Fortran Note:
3457    The Fortran interface is slightly different from that given below; it
3458    requires one to pass in  as submat a Mat (integer) array of size at least m.
3459 
3460    Level: advanced
3461 
3462    Concepts: matrices^accessing submatrices
3463    Concepts: submatrices
3464 
3465 .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal()
3466 @*/
3467 int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatReuse scall,Mat **submat)
3468 {
3469   int ierr;
3470 
3471   PetscFunctionBegin;
3472   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3473   PetscValidType(mat);
3474   MatPreallocated(mat);
3475   if (!mat->ops->getsubmatrices) SETERRQ(PETSC_ERR_SUP,"");
3476   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3477 
3478   ierr = PLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);CHKERRQ(ierr);
3479   ierr = (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);CHKERRQ(ierr);
3480   ierr = PLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);CHKERRQ(ierr);
3481 
3482   PetscFunctionReturn(0);
3483 }
3484 
3485 #undef __FUNC__
3486 #define __FUNC__ /*<a name=""></a>*/"MatDestroyMatrices"
3487 /*@C
3488    MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices().
3489 
3490    Collective on Mat
3491 
3492    Input Parameters:
3493 +  n - the number of local matrices
3494 -  mat - the matrices
3495 
3496    Level: advanced
3497 
3498 .seealso: MatGetSubMatrices()
3499 @*/
3500 int MatDestroyMatrices(int n,Mat **mat)
3501 {
3502   int ierr,i;
3503 
3504   PetscFunctionBegin;
3505   if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of matrices %d",n);
3506   PetscValidPointer(mat);
3507   for (i=0; i<n; i++) {
3508     ierr = MatDestroy((*mat)[i]);CHKERRQ(ierr);
3509   }
3510   /* memory is allocated even if n = 0 */
3511   ierr = PetscFree(*mat);CHKERRQ(ierr);
3512   PetscFunctionReturn(0);
3513 }
3514 
3515 #undef __FUNC__
3516 #define __FUNC__ /*<a name=""></a>*/"MatIncreaseOverlap"
3517 /*@
3518    MatIncreaseOverlap - Given a set of submatrices indicated by index sets,
3519    replaces the index sets by larger ones that represent submatrices with
3520    additional overlap.
3521 
3522    Collective on Mat
3523 
3524    Input Parameters:
3525 +  mat - the matrix
3526 .  n   - the number of index sets
3527 .  is  - the array of pointers to index sets
3528 -  ov  - the additional overlap requested
3529 
3530    Level: developer
3531 
3532    Concepts: overlap
3533    Concepts: ASM^computing overlap
3534 
3535 .seealso: MatGetSubMatrices()
3536 @*/
3537 int MatIncreaseOverlap(Mat mat,int n,IS *is,int ov)
3538 {
3539   int ierr;
3540 
3541   PetscFunctionBegin;
3542   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3543   PetscValidType(mat);
3544   MatPreallocated(mat);
3545   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3546   if (mat->factor)     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3547 
3548   if (!ov) PetscFunctionReturn(0);
3549   if (!mat->ops->increaseoverlap) SETERRQ(PETSC_ERR_SUP,"");
3550   ierr = PLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);CHKERRQ(ierr);
3551   ierr = (*mat->ops->increaseoverlap)(mat,n,is,ov);CHKERRQ(ierr);
3552   ierr = PLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);CHKERRQ(ierr);
3553   PetscFunctionReturn(0);
3554 }
3555 
3556 #undef __FUNC__
3557 #define __FUNC__ /*<a name=""></a>*/"MatPrintHelp"
3558 /*@
3559    MatPrintHelp - Prints all the options for the matrix.
3560 
3561    Collective on Mat
3562 
3563    Input Parameter:
3564 .  mat - the matrix
3565 
3566    Options Database Keys:
3567 +  -help - Prints matrix options
3568 -  -h - Prints matrix options
3569 
3570    Level: developer
3571 
3572 .seealso: MatCreate(), MatCreateXXX()
3573 @*/
3574 int MatPrintHelp(Mat mat)
3575 {
3576   static PetscTruth called = PETSC_FALSE;
3577   int               ierr;
3578   MPI_Comm          comm;
3579 
3580   PetscFunctionBegin;
3581   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3582   PetscValidType(mat);
3583   MatPreallocated(mat);
3584 
3585   comm = mat->comm;
3586   if (!called) {
3587     ierr = (*PetscHelpPrintf)(comm,"General matrix options:\n");CHKERRQ(ierr);
3588     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_info: view basic matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr);
3589     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_info_detailed: view detailed matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr);
3590     ierr = (*PetscHelpPrintf)(comm,"  -mat_view_draw: draw nonzero matrix structure during MatAssemblyEnd()\n");CHKERRQ(ierr);
3591     ierr = (*PetscHelpPrintf)(comm,"      -draw_pause <sec>: set seconds of display pause\n");CHKERRQ(ierr);
3592     ierr = (*PetscHelpPrintf)(comm,"      -display <name>: set alternate display\n");CHKERRQ(ierr);
3593     called = PETSC_TRUE;
3594   }
3595   if (mat->ops->printhelp) {
3596     ierr = (*mat->ops->printhelp)(mat);CHKERRQ(ierr);
3597   }
3598   PetscFunctionReturn(0);
3599 }
3600 
3601 #undef __FUNC__
3602 #define __FUNC__ /*<a name=""></a>*/"MatGetBlockSize"
3603 /*@
3604    MatGetBlockSize - Returns the matrix block size; useful especially for the
3605    block row and block diagonal formats.
3606 
3607    Not Collective
3608 
3609    Input Parameter:
3610 .  mat - the matrix
3611 
3612    Output Parameter:
3613 .  bs - block size
3614 
3615    Notes:
3616    Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG.
3617    Block row formats are MATSEQBAIJ, MATMPIBAIJ
3618 
3619    Level: intermediate
3620 
3621    Concepts: matrices^block size
3622 
3623 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag()
3624 @*/
3625 int MatGetBlockSize(Mat mat,int *bs)
3626 {
3627   int ierr;
3628 
3629   PetscFunctionBegin;
3630   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3631   PetscValidType(mat);
3632   MatPreallocated(mat);
3633   PetscValidIntPointer(bs);
3634   if (!mat->ops->getblocksize) SETERRQ(PETSC_ERR_SUP,"");
3635   ierr = (*mat->ops->getblocksize)(mat,bs);CHKERRQ(ierr);
3636   PetscFunctionReturn(0);
3637 }
3638 
3639 #undef __FUNC__
3640 #define __FUNC__ /*<a name=""></a>*/"MatGetRowIJ"
3641 /*@C
3642     MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices.
3643 
3644    Collective on Mat
3645 
3646     Input Parameters:
3647 +   mat - the matrix
3648 .   shift -  0 or 1 indicating we want the indices starting at 0 or 1
3649 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3650                 symmetrized
3651 
3652     Output Parameters:
3653 +   n - number of rows in the (possibly compressed) matrix
3654 .   ia - the row pointers
3655 .   ja - the column indices
3656 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3657 
3658     Level: developer
3659 
3660 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3661 @*/
3662 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3663 {
3664   int ierr;
3665 
3666   PetscFunctionBegin;
3667   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3668   PetscValidType(mat);
3669   MatPreallocated(mat);
3670   if (ia) PetscValidIntPointer(ia);
3671   if (ja) PetscValidIntPointer(ja);
3672   PetscValidIntPointer(done);
3673   if (!mat->ops->getrowij) *done = PETSC_FALSE;
3674   else {
3675     *done = PETSC_TRUE;
3676     ierr  = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3677   }
3678   PetscFunctionReturn(0);
3679 }
3680 
3681 #undef __FUNC__
3682 #define __FUNC__ /*<a name=""></a>*/"MatGetColumnIJ"
3683 /*@C
3684     MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices.
3685 
3686     Collective on Mat
3687 
3688     Input Parameters:
3689 +   mat - the matrix
3690 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3691 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3692                 symmetrized
3693 
3694     Output Parameters:
3695 +   n - number of columns in the (possibly compressed) matrix
3696 .   ia - the column pointers
3697 .   ja - the row indices
3698 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3699 
3700     Level: developer
3701 
3702 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3703 @*/
3704 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3705 {
3706   int ierr;
3707 
3708   PetscFunctionBegin;
3709   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3710   PetscValidType(mat);
3711   MatPreallocated(mat);
3712   if (ia) PetscValidIntPointer(ia);
3713   if (ja) PetscValidIntPointer(ja);
3714   PetscValidIntPointer(done);
3715 
3716   if (!mat->ops->getcolumnij) *done = PETSC_FALSE;
3717   else {
3718     *done = PETSC_TRUE;
3719     ierr  = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3720   }
3721   PetscFunctionReturn(0);
3722 }
3723 
3724 #undef __FUNC__
3725 #define __FUNC__ /*<a name=""></a>*/"MatRestoreRowIJ"
3726 /*@C
3727     MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with
3728     MatGetRowIJ().
3729 
3730     Collective on Mat
3731 
3732     Input Parameters:
3733 +   mat - the matrix
3734 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3735 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3736                 symmetrized
3737 
3738     Output Parameters:
3739 +   n - size of (possibly compressed) matrix
3740 .   ia - the row pointers
3741 .   ja - the column indices
3742 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3743 
3744     Level: developer
3745 
3746 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3747 @*/
3748 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3749 {
3750   int ierr;
3751 
3752   PetscFunctionBegin;
3753   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3754   PetscValidType(mat);
3755   MatPreallocated(mat);
3756   if (ia) PetscValidIntPointer(ia);
3757   if (ja) PetscValidIntPointer(ja);
3758   PetscValidIntPointer(done);
3759 
3760   if (!mat->ops->restorerowij) *done = PETSC_FALSE;
3761   else {
3762     *done = PETSC_TRUE;
3763     ierr  = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3764   }
3765   PetscFunctionReturn(0);
3766 }
3767 
3768 #undef __FUNC__
3769 #define __FUNC__ /*<a name=""></a>*/"MatRestoreColumnIJ"
3770 /*@C
3771     MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with
3772     MatGetColumnIJ().
3773 
3774     Collective on Mat
3775 
3776     Input Parameters:
3777 +   mat - the matrix
3778 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3779 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3780                 symmetrized
3781 
3782     Output Parameters:
3783 +   n - size of (possibly compressed) matrix
3784 .   ia - the column pointers
3785 .   ja - the row indices
3786 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3787 
3788     Level: developer
3789 
3790 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3791 @*/
3792 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3793 {
3794   int ierr;
3795 
3796   PetscFunctionBegin;
3797   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3798   PetscValidType(mat);
3799   MatPreallocated(mat);
3800   if (ia) PetscValidIntPointer(ia);
3801   if (ja) PetscValidIntPointer(ja);
3802   PetscValidIntPointer(done);
3803 
3804   if (!mat->ops->restorecolumnij) *done = PETSC_FALSE;
3805   else {
3806     *done = PETSC_TRUE;
3807     ierr  = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr);
3808   }
3809   PetscFunctionReturn(0);
3810 }
3811 
3812 #undef __FUNC__
3813 #define __FUNC__ /*<a name=""></a>*/"MatColoringPatch"
3814 /*@C
3815     MatColoringPatch -Used inside matrix coloring routines that
3816     use MatGetRowIJ() and/or MatGetColumnIJ().
3817 
3818     Collective on Mat
3819 
3820     Input Parameters:
3821 +   mat - the matrix
3822 .   n   - number of colors
3823 -   colorarray - array indicating color for each column
3824 
3825     Output Parameters:
3826 .   iscoloring - coloring generated using colorarray information
3827 
3828     Level: developer
3829 
3830 .seealso: MatGetRowIJ(), MatGetColumnIJ()
3831 
3832 @*/
3833 int MatColoringPatch(Mat mat,int n,int *colorarray,ISColoring *iscoloring)
3834 {
3835   int ierr;
3836 
3837   PetscFunctionBegin;
3838   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3839   PetscValidType(mat);
3840   MatPreallocated(mat);
3841   PetscValidIntPointer(colorarray);
3842 
3843   if (!mat->ops->coloringpatch) {SETERRQ(PETSC_ERR_SUP,"");}
3844   else {
3845     ierr  = (*mat->ops->coloringpatch)(mat,n,colorarray,iscoloring);CHKERRQ(ierr);
3846   }
3847   PetscFunctionReturn(0);
3848 }
3849 
3850 
3851 #undef __FUNC__
3852 #define __FUNC__ /*<a name=""></a>*/"MatSetUnfactored"
3853 /*@
3854    MatSetUnfactored - Resets a factored matrix to be treated as unfactored.
3855 
3856    Collective on Mat
3857 
3858    Input Parameter:
3859 .  mat - the factored matrix to be reset
3860 
3861    Notes:
3862    This routine should be used only with factored matrices formed by in-place
3863    factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE
3864    format).  This option can save memory, for example, when solving nonlinear
3865    systems with a matrix-free Newton-Krylov method and a matrix-based, in-place
3866    ILU(0) preconditioner.
3867 
3868    Note that one can specify in-place ILU(0) factorization by calling
3869 .vb
3870      PCType(pc,PCILU);
3871      PCILUSeUseInPlace(pc);
3872 .ve
3873    or by using the options -pc_type ilu -pc_ilu_in_place
3874 
3875    In-place factorization ILU(0) can also be used as a local
3876    solver for the blocks within the block Jacobi or additive Schwarz
3877    methods (runtime option: -sub_pc_ilu_in_place).  See the discussion
3878    of these preconditioners in the users manual for details on setting
3879    local solver options.
3880 
3881    Most users should employ the simplified SLES interface for linear solvers
3882    instead of working directly with matrix algebra routines such as this.
3883    See, e.g., SLESCreate().
3884 
3885    Level: developer
3886 
3887 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace()
3888 
3889    Concepts: matrices^unfactored
3890 
3891 @*/
3892 int MatSetUnfactored(Mat mat)
3893 {
3894   int ierr;
3895 
3896   PetscFunctionBegin;
3897   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3898   PetscValidType(mat);
3899   MatPreallocated(mat);
3900   mat->factor = 0;
3901   if (!mat->ops->setunfactored) PetscFunctionReturn(0);
3902   ierr = (*mat->ops->setunfactored)(mat);CHKERRQ(ierr);
3903   PetscFunctionReturn(0);
3904 }
3905 
3906 /*MC
3907     MatGetArrayF90 - Accesses a matrix array from Fortran90.
3908 
3909     Synopsis:
3910     MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3911 
3912     Not collective
3913 
3914     Input Parameter:
3915 .   x - matrix
3916 
3917     Output Parameters:
3918 +   xx_v - the Fortran90 pointer to the array
3919 -   ierr - error code
3920 
3921     Example of Usage:
3922 .vb
3923       Scalar, pointer xx_v(:)
3924       ....
3925       call MatGetArrayF90(x,xx_v,ierr)
3926       a = xx_v(3)
3927       call MatRestoreArrayF90(x,xx_v,ierr)
3928 .ve
3929 
3930     Notes:
3931     Not yet supported for all F90 compilers
3932 
3933     Level: advanced
3934 
3935 .seealso:  MatRestoreArrayF90(), MatGetArray(), MatRestoreArray()
3936 
3937     Concepts: matrices^accessing array
3938 
3939 M*/
3940 
3941 /*MC
3942     MatRestoreArrayF90 - Restores a matrix array that has been
3943     accessed with MatGetArrayF90().
3944 
3945     Synopsis:
3946     MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3947 
3948     Not collective
3949 
3950     Input Parameters:
3951 +   x - matrix
3952 -   xx_v - the Fortran90 pointer to the array
3953 
3954     Output Parameter:
3955 .   ierr - error code
3956 
3957     Example of Usage:
3958 .vb
3959        Scalar, pointer xx_v(:)
3960        ....
3961        call MatGetArrayF90(x,xx_v,ierr)
3962        a = xx_v(3)
3963        call MatRestoreArrayF90(x,xx_v,ierr)
3964 .ve
3965 
3966     Notes:
3967     Not yet supported for all F90 compilers
3968 
3969     Level: advanced
3970 
3971 .seealso:  MatGetArrayF90(), MatGetArray(), MatRestoreArray()
3972 
3973 M*/
3974 
3975 
3976 #undef __FUNC__
3977 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrix"
3978 /*@
3979     MatGetSubMatrix - Gets a single submatrix on the same number of processors
3980                       as the original matrix.
3981 
3982     Collective on Mat
3983 
3984     Input Parameters:
3985 +   mat - the original matrix
3986 .   isrow - rows this processor should obtain
3987 .   iscol - columns for all processors you wish to keep
3988 .   csize - number of columns "local" to this processor (does nothing for sequential
3989             matrices). This should match the result from VecGetLocalSize(x,...) if you
3990             plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE
3991 -   cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3992 
3993     Output Parameter:
3994 .   newmat - the new submatrix, of the same type as the old
3995 
3996     Level: advanced
3997 
3998     Notes: the iscol argument MST be the same on each processor
3999 
4000     Concepts: matrices^submatrices
4001 
4002 .seealso: MatGetSubMatrices()
4003 @*/
4004 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat)
4005 {
4006   int     ierr, size;
4007   Mat     *local;
4008 
4009   PetscFunctionBegin;
4010   PetscValidType(mat);
4011   MatPreallocated(mat);
4012   ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr);
4013 
4014   /* if original matrix is on just one processor then use submatrix generated */
4015   if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) {
4016     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr);
4017     PetscFunctionReturn(0);
4018   } else if (!mat->ops->getsubmatrix && size == 1) {
4019     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr);
4020     *newmat = *local;
4021     ierr = PetscFree(local);CHKERRQ(ierr);
4022     PetscFunctionReturn(0);
4023   }
4024 
4025   if (!mat->ops->getsubmatrix) SETERRQ(PETSC_ERR_SUP,"Not currently implemented");
4026   ierr = (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr);
4027   PetscFunctionReturn(0);
4028 }
4029 
4030 #undef __FUNC__
4031 #define __FUNC__ /*<a name=""></a>*/"MatGetMaps"
4032 /*@C
4033    MatGetMaps - Returns the maps associated with the matrix.
4034 
4035    Not Collective
4036 
4037    Input Parameter:
4038 .  mat - the matrix
4039 
4040    Output Parameters:
4041 +  rmap - the row (right) map
4042 -  cmap - the column (left) map
4043 
4044    Level: developer
4045 
4046    Concepts: maps^getting from matrix
4047 
4048 @*/
4049 int MatGetMaps(Mat mat,Map *rmap,Map *cmap)
4050 {
4051   int ierr;
4052 
4053   PetscFunctionBegin;
4054   PetscValidHeaderSpecific(mat,MAT_COOKIE);
4055   PetscValidType(mat);
4056   MatPreallocated(mat);
4057   ierr = (*mat->ops->getmaps)(mat,rmap,cmap);CHKERRQ(ierr);
4058   PetscFunctionReturn(0);
4059 }
4060 
4061 /*
4062       Version that works for all PETSc matrices
4063 */
4064 #undef __FUNC__
4065 #define __FUNC__ /*<a name=""></a>*/"MatGetMaps_Petsc"
4066 int MatGetMaps_Petsc(Mat mat,Map *rmap,Map *cmap)
4067 {
4068   PetscFunctionBegin;
4069   if (rmap) *rmap = mat->rmap;
4070   if (cmap) *cmap = mat->cmap;
4071   PetscFunctionReturn(0);
4072 }
4073 
4074 #undef __FUNC__
4075 #define __FUNC__ /*<a name=""></a>*/"MatSetStashInitialSize"
4076 /*@
4077    MatSetStashInitialSize - sets the sizes of the matrix stash, that is
4078    used during the assembly process to store values that belong to
4079    other processors.
4080 
4081    Not Collective
4082 
4083    Input Parameters:
4084 +  mat   - the matrix
4085 .  size  - the initial size of the stash.
4086 -  bsize - the initial size of the block-stash(if used).
4087 
4088    Options Database Keys:
4089 +   -matstash_initial_size <size> or <size0,size1,...sizep-1>
4090 -   -matstash_block_initial_size <bsize>  or <bsize0,bsize1,...bsizep-1>
4091 
4092    Level: intermediate
4093 
4094    Notes:
4095      The block-stash is used for values set with VecSetValuesBlocked() while
4096      the stash is used for values set with VecSetValues()
4097 
4098      Run with the option -log_info and look for output of the form
4099      MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
4100      to determine the appropriate value, MM, to use for size and
4101      MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
4102      to determine the value, BMM to use for bsize
4103 
4104    Concepts: stash^setting matrix size
4105    Concepts: matrices^stash
4106 
4107 @*/
4108 int MatSetStashInitialSize(Mat mat,int size, int bsize)
4109 {
4110   int ierr;
4111 
4112   PetscFunctionBegin;
4113   PetscValidHeaderSpecific(mat,MAT_COOKIE);
4114   PetscValidType(mat);
4115   MatPreallocated(mat);
4116   ierr = MatStashSetInitialSize_Private(&mat->stash,size);CHKERRQ(ierr);
4117   ierr = MatStashSetInitialSize_Private(&mat->bstash,bsize);CHKERRQ(ierr);
4118   PetscFunctionReturn(0);
4119 }
4120 
4121 #undef __FUNC__
4122 #define __FUNC__ /*<a name=""></a>*/"MatInterpolateAdd"
4123 /*@
4124    MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of
4125      the matrix
4126 
4127    Collective on Mat
4128 
4129    Input Parameters:
4130 +  mat   - the matrix
4131 .  x,y - the vectors
4132 -  w - where the result is stored
4133 
4134    Level: intermediate
4135 
4136    Notes:
4137     w may be the same vector as y.
4138 
4139     This allows one to use either the restriction or interpolation (its transpose)
4140     matrix to do the interpolation
4141 
4142     Concepts: interpolation
4143 
4144 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
4145 
4146 @*/
4147 int MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w)
4148 {
4149   int M,N,ierr;
4150 
4151   PetscFunctionBegin;
4152   PetscValidType(A);
4153   MatPreallocated(A);
4154   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
4155   if (N > M) {
4156     ierr = MatMultTransposeAdd(A,x,y,w);CHKERRQ(ierr);
4157   } else {
4158     ierr = MatMultAdd(A,x,y,w);CHKERRQ(ierr);
4159   }
4160   PetscFunctionReturn(0);
4161 }
4162 
4163 #undef __FUNC__
4164 #define __FUNC__ /*<a name=""></a>*/"MatInterpolate"
4165 /*@
4166    MatInterpolate - y = A*x or A'*x depending on the shape of
4167      the matrix
4168 
4169    Collective on Mat
4170 
4171    Input Parameters:
4172 +  mat   - the matrix
4173 -  x,y - the vectors
4174 
4175    Level: intermediate
4176 
4177    Notes:
4178     This allows one to use either the restriction or interpolation (its transpose)
4179     matrix to do the interpolation
4180 
4181    Concepts: matrices^interpolation
4182 
4183 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()
4184 
4185 @*/
4186 int MatInterpolate(Mat A,Vec x,Vec y)
4187 {
4188   int M,N,ierr;
4189 
4190   PetscFunctionBegin;
4191   PetscValidType(A);
4192   MatPreallocated(A);
4193   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
4194   if (N > M) {
4195     ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
4196   } else {
4197     ierr = MatMult(A,x,y);CHKERRQ(ierr);
4198   }
4199   PetscFunctionReturn(0);
4200 }
4201 
4202 #undef __FUNC__
4203 #define __FUNC__ /*<a name=""></a>*/"MatRestrict"
4204 /*@
4205    MatRestrict - y = A*x or A'*x
4206 
4207    Collective on Mat
4208 
4209    Input Parameters:
4210 +  mat   - the matrix
4211 -  x,y - the vectors
4212 
4213    Level: intermediate
4214 
4215    Notes:
4216     This allows one to use either the restriction or interpolation (its transpose)
4217     matrix to do the restriction
4218 
4219    Concepts: matrices^restriction
4220 
4221 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatInterpolate()
4222 
4223 @*/
4224 int MatRestrict(Mat A,Vec x,Vec y)
4225 {
4226   int M,N,ierr;
4227 
4228   PetscFunctionBegin;
4229   PetscValidType(A);
4230   MatPreallocated(A);
4231   ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
4232   if (N > M) {
4233     ierr = MatMult(A,x,y);CHKERRQ(ierr);
4234   } else {
4235     ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);
4236   }
4237   PetscFunctionReturn(0);
4238 }
4239 
4240 #undef __FUNC__
4241 #define __FUNC__ /*<a name=""></a>*/"MatNullSpaceAttach"
4242 /*@C
4243    MatNullSpaceAttach - attaches a null space to a matrix.
4244         This null space will be removed from the resulting vector whenever
4245         MatMult() is called
4246 
4247    Collective on Mat
4248 
4249    Input Parameters:
4250 +  mat - the matrix
4251 -  nullsp - the null space object
4252 
4253    Level: developer
4254 
4255    Notes:
4256       Overwrites any previous null space that may have been attached
4257 
4258    Concepts: null space^attaching to matrix
4259 
4260 .seealso: MatCreate(), MatNullSpaceCreate()
4261 @*/
4262 int MatNullSpaceAttach(Mat mat,MatNullSpace nullsp)
4263 {
4264   int ierr = 0;
4265 
4266   PetscFunctionBegin;
4267   PetscValidHeaderSpecific(mat,MAT_COOKIE);
4268   PetscValidType(mat);
4269   MatPreallocated(mat);
4270   PetscValidHeaderSpecific(nullsp,MATNULLSPACE_COOKIE);
4271 
4272   if (mat->nullsp) {
4273     ierr = MatNullSpaceDestroy(mat->nullsp);CHKERRQ(ierr);
4274   }
4275   mat->nullsp = nullsp;
4276   ierr = PetscObjectReference((PetscObject)nullsp);CHKERRQ(ierr);
4277   PetscFunctionReturn(0);
4278 }
4279 
4280 #undef __FUNC__
4281 #define __FUNC__ /*<a name=""></a>*/"MatIncompleteCholeskyFactor"
4282 /*@
4283    MatIncompleteCholeskyFactor - Performs in-place incomplete Cholesky factorization of matrix.
4284 
4285    Collective on Mat
4286 
4287    Input Parameters:
4288 +  mat - the matrix
4289 .  row - row/column permutation
4290 .  fill - expected fill factor >= 1.0
4291 -  level - level of fill, for ICC(k)
4292 
4293    Notes:
4294    Probably really in-place only when level of fill is zero, otherwise allocates
4295    new space to store factored matrix and deletes previous memory.
4296 
4297    Most users should employ the simplified SLES interface for linear solvers
4298    instead of working directly with matrix algebra routines such as this.
4299    See, e.g., SLESCreate().
4300 
4301    Level: developer
4302 
4303    Concepts: matrices^incomplete Cholesky factorization
4304    Concepts: Cholesky factorization
4305 
4306 .seealso: MatIncompleteCholeskyFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
4307 @*/
4308 int MatIncompleteCholeskyFactor(Mat mat,IS row,PetscReal fill,int level)
4309 {
4310   int ierr;
4311 
4312   PetscFunctionBegin;
4313   PetscValidHeaderSpecific(mat,MAT_COOKIE);
4314   PetscValidType(mat);
4315   MatPreallocated(mat);
4316   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
4317   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4318   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4319   if (!mat->ops->incompletecholeskyfactor) SETERRQ(PETSC_ERR_SUP,"");
4320 
4321   ierr = (*mat->ops->incompletecholeskyfactor)(mat,row,fill,level);CHKERRQ(ierr);
4322   PetscFunctionReturn(0);
4323 }
4324