xref: /petsc/src/mat/interface/matrix.c (revision 8692ea41458b8562d1c255739c878b83206e9002)
1 #ifdef PETSC_RCS_HEADER
2 static char vcid[] = "$Id: matrix.c,v 1.318 1999/02/15 21:56:30 balay Exp balay $";
3 #endif
4 
5 /*
6    This is where the abstract matrix operations are defined
7 */
8 
9 #include "src/mat/matimpl.h"        /*I "mat.h" I*/
10 #include "src/vec/vecimpl.h"
11 
12 #undef __FUNC__
13 #define __FUNC__ "MatGetRow"
14 /*@C
15    MatGetRow - Gets a row of a matrix.  You MUST call MatRestoreRow()
16    for each row that you get to ensure that your application does
17    not bleed memory.
18 
19    Not Collective
20 
21    Input Parameters:
22 +  mat - the matrix
23 -  row - the row to get
24 
25    Output Parameters:
26 +  ncols -  the number of nonzeros in the row
27 .  cols - if nonzero, the column numbers
28 -  vals - if nonzero, the values
29 
30    Notes:
31    This routine is provided for people who need to have direct access
32    to the structure of a matrix.  We hope that we provide enough
33    high-level matrix routines that few users will need it.
34 
35    MatGetRow() always returns 0-based column indices, regardless of
36    whether the internal representation is 0-based (default) or 1-based.
37 
38    For better efficiency, set cols and/or vals to PETSC_NULL if you do
39    not wish to extract these quantities.
40 
41    The user can only examine the values extracted with MatGetRow();
42    the values cannot be altered.  To change the matrix entries, one
43    must use MatSetValues().
44 
45    You can only have one call to MatGetRow() outstanding for a particular
46    matrix at a time.
47 
48    Fortran Notes:
49    The calling sequence from Fortran is
50 .vb
51    MatGetRow(matrix,row,ncols,cols,values,ierr)
52          Mat     matrix (input)
53          integer row    (input)
54          integer ncols  (output)
55          integer cols(maxcols) (output)
56          double precision (or double complex) values(maxcols) output
57 .ve
58    where maxcols >= maximum nonzeros in any row of the matrix.
59 
60    Caution:
61    Do not try to change the contents of the output arrays (cols and vals).
62    In some cases, this may corrupt the matrix.
63 
64 .keywords: matrix, row, get, extract
65 
66 .seealso: MatRestoreRow(), MatSetValues()
67 @*/
68 int MatGetRow(Mat mat,int row,int *ncols,int **cols,Scalar **vals)
69 {
70   int   ierr;
71 
72   PetscFunctionBegin;
73   PetscValidHeaderSpecific(mat,MAT_COOKIE);
74   PetscValidIntPointer(ncols);
75   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
76   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
77   if (!mat->ops->getrow) SETERRQ(PETSC_ERR_SUP,0,"");
78   PLogEventBegin(MAT_GetRow,mat,0,0,0);
79   ierr = (*mat->ops->getrow)(mat,row,ncols,cols,vals); CHKERRQ(ierr);
80   PLogEventEnd(MAT_GetRow,mat,0,0,0);
81   PetscFunctionReturn(0);
82 }
83 
84 #undef __FUNC__
85 #define __FUNC__ "MatRestoreRow"
86 /*@C
87    MatRestoreRow - Frees any temporary space allocated by MatGetRow().
88 
89    Not Collective
90 
91    Input Parameters:
92 +  mat - the matrix
93 .  row - the row to get
94 .  ncols, cols - the number of nonzeros and their columns
95 -  vals - if nonzero the column values
96 
97    Notes:
98    This routine should be called after you have finished examining the entries.
99 
100    Fortran Notes:
101    The calling sequence from Fortran is
102 .vb
103    MatRestoreRow(matrix,row,ncols,cols,values,ierr)
104       Mat     matrix (input)
105       integer row    (input)
106       integer ncols  (output)
107       integer cols(maxcols) (output)
108       double precision (or double complex) values(maxcols) output
109 .ve
110    Where maxcols >= maximum nonzeros in any row of the matrix.
111 
112    In Fortran MatRestoreRow() MUST be called after MatGetRow()
113    before another call to MatGetRow() can be made.
114 
115 .keywords: matrix, row, restore
116 
117 .seealso:  MatGetRow()
118 @*/
119 int MatRestoreRow(Mat mat,int row,int *ncols,int **cols,Scalar **vals)
120 {
121   int ierr;
122 
123   PetscFunctionBegin;
124   PetscValidHeaderSpecific(mat,MAT_COOKIE);
125   PetscValidIntPointer(ncols);
126   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
127   if (!mat->ops->restorerow) PetscFunctionReturn(0);
128   ierr = (*mat->ops->restorerow)(mat,row,ncols,cols,vals);CHKERRQ(ierr);
129   PetscFunctionReturn(0);
130 }
131 
132 #undef __FUNC__
133 #define __FUNC__ "MatView"
134 /*@C
135    MatView - Visualizes a matrix object.
136 
137    Collective on Mat unless Viewer is VIEWER_STDOUT_SELF
138 
139    Input Parameters:
140 +  mat - the matrix
141 -  ptr - visualization context
142 
143   Notes:
144   The available visualization contexts include
145 +    VIEWER_STDOUT_SELF - standard output (default)
146 .    VIEWER_STDOUT_WORLD - synchronized standard
147         output where only the first processor opens
148         the file.  All other processors send their
149         data to the first processor to print.
150 -     VIEWER_DRAW_WORLD - graphical display of nonzero structure
151 
152    The user can open alternative visualization contexts with
153 +    ViewerASCIIOpen() - Outputs matrix to a specified file
154 .    ViewerBinaryOpen() - Outputs matrix in binary to a
155          specified file; corresponding input uses MatLoad()
156 .    ViewerDrawOpen() - Outputs nonzero matrix structure to
157          an X window display
158 -    ViewerSocketOpen() - Outputs matrix to Socket viewer.
159          Currently only the sequential dense and AIJ
160          matrix types support the Socket viewer.
161 
162    The user can call ViewerSetFormat() to specify the output
163    format of ASCII printed objects (when using VIEWER_STDOUT_SELF,
164    VIEWER_STDOUT_WORLD and ViewerASCIIOpen).  Available formats include
165 +    VIEWER_FORMAT_ASCII_DEFAULT - default, prints matrix contents
166 .    VIEWER_FORMAT_ASCII_MATLAB - prints matrix contents in Matlab format
167 .    VIEWER_FORMAT_ASCII_DENSE - prints entire matrix including zeros
168 .    VIEWER_FORMAT_ASCII_COMMON - prints matrix contents, using a sparse
169          format common among all matrix types
170 .    VIEWER_FORMAT_ASCII_IMPL - prints matrix contents, using an implementation-specific
171          format (which is in many cases the same as the default)
172 .    VIEWER_FORMAT_ASCII_INFO - prints basic information about the matrix
173          size and structure (not the matrix entries)
174 -    VIEWER_FORMAT_ASCII_INFO_LONG - prints more detailed information about
175          the matrix structure
176 
177 .keywords: matrix, view, visualize, output, print, write, draw
178 
179 .seealso: ViewerSetFormat(), ViewerASCIIOpen(), ViewerDrawOpen(),
180           ViewerSocketOpen(), ViewerBinaryOpen(), MatLoad()
181 @*/
182 int MatView(Mat mat,Viewer viewer)
183 {
184   int          format, ierr, rows, cols;
185   char         *cstr;
186   ViewerType   vtype;
187 
188   PetscFunctionBegin;
189   PetscValidHeaderSpecific(mat,MAT_COOKIE);
190   if (viewer) PetscValidHeaderSpecific(viewer,VIEWER_COOKIE);
191   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
192 
193   if (!viewer) {
194     viewer = VIEWER_STDOUT_SELF;
195   }
196 
197   ierr = ViewerGetType(viewer,&vtype);
198   if (PetscTypeCompare(vtype,ASCII_VIEWER)) {
199     ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr);
200     if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) {
201       ViewerASCIIPrintf(viewer,"Matrix Object:\n");
202       ierr = MatGetType(mat,PETSC_NULL,&cstr); CHKERRQ(ierr);
203       ierr = MatGetSize(mat,&rows,&cols); CHKERRQ(ierr);
204       ViewerASCIIPrintf(viewer,"  type=%s, rows=%d, cols=%d\n",cstr,rows,cols);
205       if (mat->ops->getinfo) {
206         MatInfo info;
207         ierr = MatGetInfo(mat,MAT_GLOBAL_SUM,&info); CHKERRQ(ierr);
208         ViewerASCIIPrintf(viewer,"  total: nonzeros=%d, allocated nonzeros=%d\n",
209                           (int)info.nz_used,(int)info.nz_allocated);
210       }
211     }
212   }
213   if (mat->ops->view) {
214     ierr = ViewerASCIIPushTab(viewer);CHKERRQ(ierr);
215     ierr = (*mat->ops->view)(mat,viewer); CHKERRQ(ierr);
216     ierr = ViewerASCIIPopTab(viewer);CHKERRQ(ierr);
217   }
218   PetscFunctionReturn(0);
219 }
220 
221 #undef __FUNC__
222 #define __FUNC__ "MatScaleSystem"
223 /*@C
224    MatScaleSystem - Scale a vector solution and right hand side to
225      match the scaling of a scaled matrix.
226 
227    Collective on Mat
228 
229    Input Parameter:
230 +  mat - the matrix
231 .  x - solution vector (or PETSC_NULL)
232 +  b - right hand side vector (or PETSC_NULL)
233 
234 
235   Notes:
236     For AIJ, BAIJ, and BDiag matrices the matrices are not
237   internally scaled, so this does nothing. For MPIRowBS it
238   permutes and diagonally scales.
239 
240     The SLES methods automatically call this routine when required
241   (via PCPreSolve()) so it is rarely used directly.
242 
243   Level: Developer
244 
245 .keywords: matrix, scale
246 
247 .seealso: MatUseScaledForm(), MatUnScaleSystem()
248 @*/
249 int MatScaleSystem(Mat mat,Vec x,Vec b)
250 {
251   int ierr;
252 
253   PetscFunctionBegin;
254   PetscValidHeaderSpecific(mat,MAT_COOKIE);
255   if (x) {PetscValidHeaderSpecific(x,VEC_COOKIE);}
256   if (b) {PetscValidHeaderSpecific(b,VEC_COOKIE);}
257   if (mat->ops->scalesystem) {
258     ierr = (*mat->ops->scalesystem)(mat,x,b); CHKERRQ(ierr);
259   }
260   PetscFunctionReturn(0);
261 }
262 
263 #undef __FUNC__
264 #define __FUNC__ "MatUnScaleSystem"
265 /*@C
266    MatUnScaleSystem - Unscales a vector solution and right hand side to
267      match the original scaling of a scaled matrix.
268 
269    Collective on Mat
270 
271    Input Parameter:
272 +  mat - the matrix
273 .  x - solution vector (or PETSC_NULL)
274 +  b - right hand side vector (or PETSC_NULL)
275 
276 
277   Notes:
278     For AIJ, BAIJ, and BDiag matrices the matrices are not
279   internally scaled, so this does nothing. For MPIRowBS it
280   permutes and diagonally scales.
281 
282     The SLES methods automatically call this routine when required
283   (via PCPostSolve()) so it is rarely used directly.
284 
285   Level: Developer
286 
287 .keywords: matrix, scale
288 
289 .seealso: MatUseScaledForm(), MatScaleSystem()
290 @*/
291 int MatUnScaleSystem(Mat mat,Vec x,Vec b)
292 {
293   int ierr;
294 
295   PetscFunctionBegin;
296   PetscValidHeaderSpecific(mat,MAT_COOKIE);
297   if (x) {PetscValidHeaderSpecific(x,VEC_COOKIE);}
298   if (b) {PetscValidHeaderSpecific(b,VEC_COOKIE);}
299   if (mat->ops->unscalesystem) {
300     ierr = (*mat->ops->unscalesystem)(mat,x,b); CHKERRQ(ierr);
301   }
302   PetscFunctionReturn(0);
303 }
304 
305 #undef __FUNC__
306 #define __FUNC__ "MatUseScaledForm"
307 /*@C
308    MatUseScaledForm - For matrix storage formats that scale the
309      matrix (for example MPIRowbs matrices are diagonally scaled on
310      assembly) indicates matrix operations (MatMult() etc) are
311      applied using the scaled matrix.
312 
313    Collective on Mat
314 
315    Input Parameter:
316 +  mat - the matrix
317 -  scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for
318             applying the original matrix
319 
320   Notes:
321     For scaled formats applying the original, unscaled matrix
322    will be slightly more expensive
323 
324   Level: Developer
325 
326 .keywords: matrix, scale
327 
328 .seealso: MatScaleSystem(), MatUnScaleSystem()
329 @*/
330 int MatUseScaledForm(Mat mat,PetscTruth scaled)
331 {
332   int ierr;
333 
334   PetscFunctionBegin;
335   PetscValidHeaderSpecific(mat,MAT_COOKIE);
336   if (mat->ops->usescaledform) {
337     ierr = (*mat->ops->usescaledform)(mat,scaled); CHKERRQ(ierr);
338   }
339   PetscFunctionReturn(0);
340 }
341 
342 #undef __FUNC__
343 #define __FUNC__ "MatDestroy"
344 /*@C
345    MatDestroy - Frees space taken by a matrix.
346 
347    Collective on Mat
348 
349    Input Parameter:
350 .  mat - the matrix
351 
352 .keywords: matrix, destroy
353 @*/
354 int MatDestroy(Mat mat)
355 {
356   int ierr;
357 
358   PetscFunctionBegin;
359   PetscValidHeaderSpecific(mat,MAT_COOKIE);
360   ierr = (*mat->ops->destroy)(mat); CHKERRQ(ierr);
361   PetscFunctionReturn(0);
362 }
363 
364 #undef __FUNC__
365 #define __FUNC__ "MatValid"
366 /*@
367    MatValid - Checks whether a matrix object is valid.
368 
369    Collective on Mat
370 
371    Input Parameter:
372 .  m - the matrix to check
373 
374    Output Parameter:
375    flg - flag indicating matrix status, either
376    PETSC_TRUE if matrix is valid, or PETSC_FALSE otherwise.
377 
378 .keywords: matrix, valid
379 @*/
380 int MatValid(Mat m,PetscTruth *flg)
381 {
382   PetscFunctionBegin;
383   PetscValidIntPointer(flg);
384   if (!m)                           *flg = PETSC_FALSE;
385   else if (m->cookie != MAT_COOKIE) *flg = PETSC_FALSE;
386   else                              *flg = PETSC_TRUE;
387   PetscFunctionReturn(0);
388 }
389 
390 #undef __FUNC__
391 #define __FUNC__ "MatSetValues"
392 /*@
393    MatSetValues - Inserts or adds a block of values into a matrix.
394    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
395    MUST be called after all calls to MatSetValues() have been completed.
396 
397    Not Collective
398 
399    Input Parameters:
400 +  mat - the matrix
401 .  v - a logically two-dimensional array of values
402 .  m, idxm - the number of rows and their global indices
403 .  n, idxn - the number of columns and their global indices
404 -  addv - either ADD_VALUES or INSERT_VALUES, where
405    ADD_VALUES adds values to any existing entries, and
406    INSERT_VALUES replaces existing entries with new values
407 
408    Notes:
409    By default the values, v, are row-oriented and unsorted.
410    See MatSetOption() for other options.
411 
412    Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES
413    options cannot be mixed without intervening calls to the assembly
414    routines.
415 
416    MatSetValues() uses 0-based row and column numbers in Fortran
417    as well as in C.
418 
419    Negative indices may be passed in idxm and idxn, these rows and columns are
420    simply ignored. This allows easily inserting element stiffness matrices
421    with homogeneous Dirchlet boundary conditions that you don't want represented
422    in the matrix.
423 
424    Efficiency Alert:
425    The routine MatSetValuesBlocked() may offer much better efficiency
426    for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
427 
428    Level: beginner
429 
430 .keywords: matrix, insert, add, set, values
431 
432 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
433 @*/
434 int MatSetValues(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v,InsertMode addv)
435 {
436   int ierr;
437 
438   PetscFunctionBegin;
439   if (!m || !n) PetscFunctionReturn(0); /* no values to insert */
440   PetscValidHeaderSpecific(mat,MAT_COOKIE);
441   PetscValidIntPointer(idxm);
442   PetscValidIntPointer(idxn);
443   PetscValidScalarPointer(v);
444   if (mat->insertmode == NOT_SET_VALUES) {
445     mat->insertmode = addv;
446   }
447 #if defined(USE_PETSC_BOPT_g)
448   else if (mat->insertmode != addv) {
449     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values");
450   }
451   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
452 #endif
453 
454   if (mat->assembled) {
455     mat->was_assembled = PETSC_TRUE;
456     mat->assembled     = PETSC_FALSE;
457   }
458   PLogEventBegin(MAT_SetValues,mat,0,0,0);
459   if (!mat->ops->setvalues) SETERRQ(PETSC_ERR_SUP,1,"Not supported for this matrix type");
460   ierr = (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr);
461   PLogEventEnd(MAT_SetValues,mat,0,0,0);
462   PetscFunctionReturn(0);
463 }
464 
465 #undef __FUNC__
466 #define __FUNC__ "MatSetValuesBlocked"
467 /*@
468    MatSetValuesBlocked - Inserts or adds a block of values into a matrix.
469 
470    Not Collective
471 
472    Input Parameters:
473 +  mat - the matrix
474 .  v - a logically two-dimensional array of values
475 .  m, idxm - the number of block rows and their global block indices
476 .  n, idxn - the number of block columns and their global block indices
477 -  addv - either ADD_VALUES or INSERT_VALUES, where
478    ADD_VALUES adds values to any existing entries, and
479    INSERT_VALUES replaces existing entries with new values
480 
481    Notes:
482    By default the values, v, are row-oriented and unsorted. So the layout of
483    v is the same as for MatSetValues(). See MatSetOption() for other options.
484 
485    Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
486    options cannot be mixed without intervening calls to the assembly
487    routines.
488 
489    MatSetValuesBlocked() uses 0-based row and column numbers in Fortran
490    as well as in C.
491 
492    Negative indices may be passed in idxm and idxn, these rows and columns are
493    simply ignored. This allows easily inserting element stiffness matrices
494    with homogeneous Dirchlet boundary conditions that you don't want represented
495    in the matrix.
496 
497    Each time an entry is set within a sparse matrix via MatSetValues(),
498    internal searching must be done to determine where to place the the
499    data in the matrix storage space.  By instead inserting blocks of
500    entries via MatSetValuesBlocked(), the overhead of matrix assembly is
501    reduced.
502 
503    Restrictions:
504    MatSetValuesBlocked() is currently supported only for the block AIJ
505    matrix format (MATSEQBAIJ and MATMPIBAIJ, which are created via
506    MatCreateSeqBAIJ() and MatCreateMPIBAIJ()).
507 
508 .keywords: matrix, insert, add, set, values
509 
510 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
511 @*/
512 int MatSetValuesBlocked(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v,InsertMode addv)
513 {
514   int ierr;
515 
516   PetscFunctionBegin;
517   if (!m || !n) PetscFunctionReturn(0); /* no values to insert */
518   PetscValidHeaderSpecific(mat,MAT_COOKIE);
519   PetscValidIntPointer(idxm);
520   PetscValidIntPointer(idxn);
521   PetscValidScalarPointer(v);
522   if (mat->insertmode == NOT_SET_VALUES) {
523     mat->insertmode = addv;
524   }
525 #if defined(USE_PETSC_BOPT_g)
526   else if (mat->insertmode != addv) {
527     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values");
528   }
529   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
530 #endif
531 
532   if (mat->assembled) {
533     mat->was_assembled = PETSC_TRUE;
534     mat->assembled     = PETSC_FALSE;
535   }
536   PLogEventBegin(MAT_SetValues,mat,0,0,0);
537   if (!mat->ops->setvaluesblocked) SETERRQ(PETSC_ERR_SUP,1,"Not supported for this matrix type");
538   ierr = (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr);
539   PLogEventEnd(MAT_SetValues,mat,0,0,0);
540   PetscFunctionReturn(0);
541 }
542 
543 /*MC
544    MatSetValue - Set a single entry into a matrix.
545 
546    Synopsis:
547    void MatSetValue(Mat m,int row,int col,Scalar value,InsertMode mode);
548 
549    Not collective
550 
551    Input Parameters:
552 +  m - the matrix
553 .  row - the row location of the entry
554 .  col - the column location of the entry
555 .  value - the value to insert
556 -  mode - either INSERT_VALUES or ADD_VALUES
557 
558    Notes:
559    For efficiency one should use MatSetValues() and set several or many
560    values simultaneously if possible.
561 
562    Note that VecSetValue() does NOT return an error code (since this
563    is checked internally).
564 
565 .seealso: MatSetValues()
566 M*/
567 
568 #undef __FUNC__
569 #define __FUNC__ "MatGetValues"
570 /*@
571    MatGetValues - Gets a block of values from a matrix.
572 
573    Not Collective; currently only returns a local block
574 
575    Input Parameters:
576 +  mat - the matrix
577 .  v - a logically two-dimensional array for storing the values
578 .  m, idxm - the number of rows and their global indices
579 -  n, idxn - the number of columns and their global indices
580 
581    Notes:
582    The user must allocate space (m*n Scalars) for the values, v.
583    The values, v, are then returned in a row-oriented format,
584    analogous to that used by default in MatSetValues().
585 
586    MatGetValues() uses 0-based row and column numbers in
587    Fortran as well as in C.
588 
589    MatGetValues() requires that the matrix has been assembled
590    with MatAssemblyBegin()/MatAssemblyEnd().  Thus, calls to
591    MatSetValues() and MatGetValues() CANNOT be made in succession
592    without intermediate matrix assembly.
593 
594 .keywords: matrix, get, values
595 
596 .seealso: MatGetRow(), MatGetSubmatrices(), MatSetValues()
597 @*/
598 int MatGetValues(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v)
599 {
600   int ierr;
601 
602   PetscFunctionBegin;
603   PetscValidHeaderSpecific(mat,MAT_COOKIE);
604   PetscValidIntPointer(idxm);
605   PetscValidIntPointer(idxn);
606   PetscValidScalarPointer(v);
607   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
608   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
609   if (!mat->ops->getvalues) SETERRQ(PETSC_ERR_SUP,0,"");
610 
611   PLogEventBegin(MAT_GetValues,mat,0,0,0);
612   ierr = (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v); CHKERRQ(ierr);
613   PLogEventEnd(MAT_GetValues,mat,0,0,0);
614   PetscFunctionReturn(0);
615 }
616 
617 #undef __FUNC__
618 #define __FUNC__ "MatSetLocalToGlobalMapping"
619 /*@
620    MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
621    the routine MatSetValuesLocal() to allow users to insert matrix entries
622    using a local (per-processor) numbering.
623 
624    Not Collective
625 
626    Input Parameters:
627 +  x - the matrix
628 -  mapping - mapping created with ISLocalToGlobalMappingCreate()
629              or ISLocalToGlobalMappingCreateIS()
630 
631 .keywords: matrix, set, values, local, global, mapping
632 
633 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
634 @*/
635 int MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping)
636 {
637   PetscFunctionBegin;
638   PetscValidHeaderSpecific(x,MAT_COOKIE);
639   PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE);
640 
641   if (x->mapping) {
642     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Mapping already set for matrix");
643   }
644 
645   x->mapping = mapping;
646   PetscObjectReference((PetscObject)mapping);
647   PetscFunctionReturn(0);
648 }
649 
650 #undef __FUNC__
651 #define __FUNC__ "MatSetLocalToGlobalMappingBlocked"
652 /*@
653    MatSetLocalToGlobalMappingBlocked - Sets a local-to-global numbering for use
654    by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
655    entries using a local (per-processor) numbering.
656 
657    Not Collective
658 
659    Input Parameters:
660 +  x - the matrix
661 -  mapping - mapping created with ISLocalToGlobalMappingCreate() or
662              ISLocalToGlobalMappingCreateIS()
663 
664 .keywords: matrix, set, values, local ordering
665 
666 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
667            MatSetValuesBlocked(), MatSetValuesLocal()
668 @*/
669 int MatSetLocalToGlobalMappingBlocked(Mat x,ISLocalToGlobalMapping mapping)
670 {
671   PetscFunctionBegin;
672   PetscValidHeaderSpecific(x,MAT_COOKIE);
673   PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE);
674 
675   if (x->bmapping) {
676     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Mapping already set for matrix");
677   }
678 
679   x->bmapping = mapping;
680   PetscObjectReference((PetscObject)mapping);
681   PetscFunctionReturn(0);
682 }
683 
684 #undef __FUNC__
685 #define __FUNC__ "MatSetValuesLocal"
686 /*@
687    MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
688    using a local ordering of the nodes.
689 
690    Not Collective
691 
692    Input Parameters:
693 +  x - the matrix
694 .  nrow, irow - number of rows and their local indices
695 .  ncol, icol - number of columns and their local indices
696 .  y -  a logically two-dimensional array of values
697 -  addv - either INSERT_VALUES or ADD_VALUES, where
698    ADD_VALUES adds values to any existing entries, and
699    INSERT_VALUES replaces existing entries with new values
700 
701    Notes:
702    Before calling MatSetValuesLocal(), the user must first set the
703    local-to-global mapping by calling MatSetLocalToGlobalMapping().
704 
705    Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES
706    options cannot be mixed without intervening calls to the assembly
707    routines.
708 
709    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
710    MUST be called after all calls to MatSetValuesLocal() have been completed.
711 
712 .keywords: matrix, set, values, local ordering
713 
714 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping()
715 @*/
716 int MatSetValuesLocal(Mat mat,int nrow,int *irow,int ncol, int *icol,Scalar *y,InsertMode addv)
717 {
718   int ierr,irowm[2048],icolm[2048];
719 
720   PetscFunctionBegin;
721   PetscValidHeaderSpecific(mat,MAT_COOKIE);
722   PetscValidIntPointer(irow);
723   PetscValidIntPointer(icol);
724   PetscValidScalarPointer(y);
725 
726   if (mat->insertmode == NOT_SET_VALUES) {
727     mat->insertmode = addv;
728   }
729 #if defined(USE_PETSC_BOPT_g)
730   else if (mat->insertmode != addv) {
731     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values");
732   }
733   if (!mat->mapping) {
734     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Local to global never set with MatSetLocalToGlobalMapping()");
735   }
736   if (nrow > 2048 || ncol > 2048) {
737     SETERRQ2(PETSC_ERR_SUP,0,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
738   }
739   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
740 #endif
741 
742   if (mat->assembled) {
743     mat->was_assembled = PETSC_TRUE;
744     mat->assembled     = PETSC_FALSE;
745   }
746   PLogEventBegin(MAT_SetValues,mat,0,0,0);
747   ierr = ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm); CHKERRQ(ierr);
748   ierr = ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);CHKERRQ(ierr);
749   ierr = (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr);
750   PLogEventEnd(MAT_SetValues,mat,0,0,0);
751   PetscFunctionReturn(0);
752 }
753 
754 #undef __FUNC__
755 #define __FUNC__ "MatSetValuesBlockedLocal"
756 /*@
757    MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
758    using a local ordering of the nodes a block at a time.
759 
760    Not Collective
761 
762    Input Parameters:
763 +  x - the matrix
764 .  nrow, irow - number of rows and their local indices
765 .  ncol, icol - number of columns and their local indices
766 .  y -  a logically two-dimensional array of values
767 -  addv - either INSERT_VALUES or ADD_VALUES, where
768    ADD_VALUES adds values to any existing entries, and
769    INSERT_VALUES replaces existing entries with new values
770 
771    Notes:
772    Before calling MatSetValuesBlockedLocal(), the user must first set the
773    local-to-global mapping by calling MatSetLocalToGlobalMappingBlocked(),
774    where the mapping MUST be set for matrix blocks, not for matrix elements.
775 
776    Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
777    options cannot be mixed without intervening calls to the assembly
778    routines.
779 
780    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
781    MUST be called after all calls to MatSetValuesBlockedLocal() have been completed.
782 
783 .keywords: matrix, set, values, blocked, local
784 
785 .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlocked(), MatSetValuesBlocked()
786 @*/
787 int MatSetValuesBlockedLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,Scalar *y,InsertMode addv)
788 {
789   int ierr,irowm[2048],icolm[2048];
790 
791   PetscFunctionBegin;
792   PetscValidHeaderSpecific(mat,MAT_COOKIE);
793   PetscValidIntPointer(irow);
794   PetscValidIntPointer(icol);
795   PetscValidScalarPointer(y);
796   if (mat->insertmode == NOT_SET_VALUES) {
797     mat->insertmode = addv;
798   }
799 #if defined(USE_PETSC_BOPT_g)
800   else if (mat->insertmode != addv) {
801     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values");
802   }
803   if (!mat->bmapping) {
804     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Local to global never set with MatSetLocalToGlobalMappingBlocked()");
805   }
806   if (nrow > 2048 || ncol > 2048) {
807     SETERRQ2(PETSC_ERR_SUP,0,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
808   }
809   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
810 #endif
811 
812   if (mat->assembled) {
813     mat->was_assembled = PETSC_TRUE;
814     mat->assembled     = PETSC_FALSE;
815   }
816   PLogEventBegin(MAT_SetValues,mat,0,0,0);
817   ierr = ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm); CHKERRQ(ierr);
818   ierr = ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm); CHKERRQ(ierr);
819   ierr = (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr);
820   PLogEventEnd(MAT_SetValues,mat,0,0,0);
821   PetscFunctionReturn(0);
822 }
823 
824 /* --------------------------------------------------------*/
825 #undef __FUNC__
826 #define __FUNC__ "MatMult"
827 /*@
828    MatMult - Computes the matrix-vector product, y = Ax.
829 
830    Collective on Mat and Vec
831 
832    Input Parameters:
833 +  mat - the matrix
834 -  x   - the vector to be multilplied
835 
836    Output Parameters:
837 .  y - the result
838 
839    Notes:
840    The vectors x and y cannot be the same.  I.e., one cannot
841    call MatMult(A,y,y).
842 
843 .keywords: matrix, multiply, matrix-vector product
844 
845 .seealso: MatMultTrans(), MatMultAdd(), MatMultTransAdd()
846 @*/
847 int MatMult(Mat mat,Vec x,Vec y)
848 {
849   int ierr;
850 
851   PetscFunctionBegin;
852   PetscValidHeaderSpecific(mat,MAT_COOKIE);
853   PetscValidHeaderSpecific(x,VEC_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE);
854   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
855   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
856   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"x and y must be different vectors");
857   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
858   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
859   if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: local dim %d %d",mat->m,y->n);
860 
861   PLogEventBegin(MAT_Mult,mat,x,y,0);
862   ierr = (*mat->ops->mult)(mat,x,y); CHKERRQ(ierr);
863   PLogEventEnd(MAT_Mult,mat,x,y,0);
864 
865   PetscFunctionReturn(0);
866 }
867 
868 #undef __FUNC__
869 #define __FUNC__ "MatMultTrans"
870 /*@
871    MatMultTrans - Computes matrix transpose times a vector.
872 
873    Collective on Mat and Vec
874 
875    Input Parameters:
876 +  mat - the matrix
877 -  x   - the vector to be multilplied
878 
879    Output Parameters:
880 .  y - the result
881 
882    Notes:
883    The vectors x and y cannot be the same.  I.e., one cannot
884    call MatMultTrans(A,y,y).
885 
886 .keywords: matrix, multiply, matrix-vector product, transpose
887 
888 .seealso: MatMult(), MatMultAdd(), MatMultTransAdd()
889 @*/
890 int MatMultTrans(Mat mat,Vec x,Vec y)
891 {
892   int ierr;
893 
894   PetscFunctionBegin;
895   PetscValidHeaderSpecific(mat,MAT_COOKIE);
896   PetscValidHeaderSpecific(x,VEC_COOKIE); PetscValidHeaderSpecific(y,VEC_COOKIE);
897   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
898   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
899   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"x and y must be different vectors");
900   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
901   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
902   PLogEventBegin(MAT_MultTrans,mat,x,y,0);
903   ierr = (*mat->ops->multtrans)(mat,x,y); CHKERRQ(ierr);
904   PLogEventEnd(MAT_MultTrans,mat,x,y,0);
905   PetscFunctionReturn(0);
906 }
907 
908 #undef __FUNC__
909 #define __FUNC__ "MatMultAdd"
910 /*@
911     MatMultAdd -  Computes v3 = v2 + A * v1.
912 
913     Collective on Mat and Vec
914 
915     Input Parameters:
916 +   mat - the matrix
917 -   v1, v2 - the vectors
918 
919     Output Parameters:
920 .   v3 - the result
921 
922    Notes:
923    The vectors v1 and v3 cannot be the same.  I.e., one cannot
924    call MatMultAdd(A,v1,v2,v1).
925 
926 .keywords: matrix, multiply, matrix-vector product, add
927 
928 .seealso: MatMultTrans(), MatMult(), MatMultTransAdd()
929 @*/
930 int MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
931 {
932   int ierr;
933 
934   PetscFunctionBegin;
935   PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(v1,VEC_COOKIE);
936   PetscValidHeaderSpecific(v2,VEC_COOKIE); PetscValidHeaderSpecific(v3,VEC_COOKIE);
937   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
938   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
939   if (mat->N != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v1: global dim %d %d",mat->N,v1->N);
940   if (mat->M != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: global dim %d %d",mat->M,v2->N);
941   if (mat->M != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: global dim %d %d",mat->M,v3->N);
942   if (mat->m != v3->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: local dim %d %d",mat->m,v3->n);
943   if (mat->m != v2->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: local dim %d %d",mat->m,v2->n);
944   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,0,"v1 and v3 must be different vectors");
945 
946   PLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);
947   ierr = (*mat->ops->multadd)(mat,v1,v2,v3); CHKERRQ(ierr);
948   PLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);
949   PetscFunctionReturn(0);
950 }
951 
952 #undef __FUNC__
953 #define __FUNC__ "MatMultTransAdd"
954 /*@
955    MatMultTransAdd - Computes v3 = v2 + A' * v1.
956 
957    Collective on Mat and Vec
958 
959    Input Parameters:
960 +  mat - the matrix
961 -  v1, v2 - the vectors
962 
963    Output Parameters:
964 .  v3 - the result
965 
966    Notes:
967    The vectors v1 and v3 cannot be the same.  I.e., one cannot
968    call MatMultTransAdd(A,v1,v2,v1).
969 
970 .keywords: matrix, multiply, matrix-vector product, transpose, add
971 
972 .seealso: MatMultTrans(), MatMultAdd(), MatMult()
973 @*/
974 int MatMultTransAdd(Mat mat,Vec v1,Vec v2,Vec v3)
975 {
976   int ierr;
977 
978   PetscFunctionBegin;
979   PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(v1,VEC_COOKIE);
980   PetscValidHeaderSpecific(v2,VEC_COOKIE);PetscValidHeaderSpecific(v3,VEC_COOKIE);
981   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
982   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
983   if (!mat->ops->multtransadd) SETERRQ(PETSC_ERR_SUP,0,"");
984   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,0,"v1 and v3 must be different vectors");
985   if (mat->M != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v1: global dim %d %d",mat->M,v1->N);
986   if (mat->N != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: global dim %d %d",mat->N,v2->N);
987   if (mat->N != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: global dim %d %d",mat->N,v3->N);
988 
989   PLogEventBegin(MAT_MultTransAdd,mat,v1,v2,v3);
990   ierr = (*mat->ops->multtransadd)(mat,v1,v2,v3); CHKERRQ(ierr);
991   PLogEventEnd(MAT_MultTransAdd,mat,v1,v2,v3);
992   PetscFunctionReturn(0);
993 }
994 /* ------------------------------------------------------------*/
995 #undef __FUNC__
996 #define __FUNC__ "MatGetInfo"
997 /*@C
998    MatGetInfo - Returns information about matrix storage (number of
999    nonzeros, memory, etc.).
1000 
1001    Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used
1002    as the flag
1003 
1004    Input Parameters:
1005 .  mat - the matrix
1006 
1007    Output Parameters:
1008 +  flag - flag indicating the type of parameters to be returned
1009    (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
1010    MAT_GLOBAL_SUM - sum over all processors)
1011 -  info - matrix information context
1012 
1013    Notes:
1014    The MatInfo context contains a variety of matrix data, including
1015    number of nonzeros allocated and used, number of mallocs during
1016    matrix assembly, etc.  Additional information for factored matrices
1017    is provided (such as the fill ratio, number of mallocs during
1018    factorization, etc.).  Much of this info is printed to STDOUT
1019    when using the runtime options
1020 $       -log_info -mat_view_info
1021 
1022    Example for C/C++ Users:
1023    See the file ${PETSC_DIR}/include/mat.h for a complete list of
1024    data within the MatInfo context.  For example,
1025 .vb
1026       MatInfo info;
1027       Mat     A;
1028       double  mal, nz_a, nz_u;
1029 
1030       MatGetInfo(A,MAT_LOCAL,&info);
1031       mal  = info.mallocs;
1032       nz_a = info.nz_allocated;
1033 .ve
1034 
1035    Example for Fortran Users:
1036    Fortran users should declare info as a double precision
1037    array of dimension MAT_INFO_SIZE, and then extract the parameters
1038    of interest.  See the file ${PETSC_DIR}/include/finclude/mat.h
1039    a complete list of parameter names.
1040 .vb
1041       double  precision info(MAT_INFO_SIZE)
1042       double  precision mal, nz_a
1043       Mat     A
1044       integer ierr
1045 
1046       call MatGetInfo(A,MAT_LOCAL,info,ierr)
1047       mal = info(MAT_INFO_MALLOCS)
1048       nz_a = info(MAT_INFO_NZ_ALLOCATED)
1049 .ve
1050 
1051 .keywords: matrix, get, info, storage, nonzeros, memory, fill
1052 @*/
1053 int MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
1054 {
1055   int ierr;
1056 
1057   PetscFunctionBegin;
1058   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1059   PetscValidPointer(info);
1060   if (!mat->ops->getinfo) SETERRQ(PETSC_ERR_SUP,0,"");
1061   ierr = (*mat->ops->getinfo)(mat,flag,info);CHKERRQ(ierr);
1062   PetscFunctionReturn(0);
1063 }
1064 
1065 /* ----------------------------------------------------------*/
1066 #undef __FUNC__
1067 #define __FUNC__ "MatILUDTFactor"
1068 /*@
1069    MatILUDTFactor - Performs a drop tolerance ILU factorization.
1070 
1071    Collective on Mat
1072 
1073    Input Parameters:
1074 +  mat - the matrix
1075 .  dt  - the drop tolerance
1076 .  maxnz - the maximum number of nonzeros per row allowed?
1077 .  row - row permutation
1078 -  col - column permutation
1079 
1080    Output Parameters:
1081 .  fact - the factored matrix
1082 
1083 .keywords: matrix, factor, LU, in-place
1084 
1085 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
1086 @*/
1087 int MatILUDTFactor(Mat mat,double dt,int maxnz,IS row,IS col,Mat *fact)
1088 {
1089   int ierr;
1090 
1091   PetscFunctionBegin;
1092   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1093   PetscValidPointer(fact);
1094   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1095   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1096   if (!mat->ops->iludtfactor) SETERRQ(PETSC_ERR_SUP,0,"");
1097 
1098   PLogEventBegin(MAT_ILUFactor,mat,row,col,0);
1099   ierr = (*mat->ops->iludtfactor)(mat,dt,maxnz,row,col,fact); CHKERRQ(ierr);
1100   PLogEventEnd(MAT_ILUFactor,mat,row,col,0);
1101 
1102   PetscFunctionReturn(0);
1103 }
1104 
1105 #undef __FUNC__
1106 #define __FUNC__ "MatLUFactor"
1107 /*@
1108    MatLUFactor - Performs in-place LU factorization of matrix.
1109 
1110    Collective on Mat
1111 
1112    Input Parameters:
1113 +  mat - the matrix
1114 .  row - row permutation
1115 .  col - column permutation
1116 -  f - expected fill as ratio of original fill.
1117 
1118 .keywords: matrix, factor, LU, in-place
1119 
1120 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
1121 @*/
1122 int MatLUFactor(Mat mat,IS row,IS col,double f)
1123 {
1124   int ierr;
1125 
1126   PetscFunctionBegin;
1127   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1128   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square");
1129   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1130   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1131   if (!mat->ops->lufactor) SETERRQ(PETSC_ERR_SUP,0,"");
1132 
1133   PLogEventBegin(MAT_LUFactor,mat,row,col,0);
1134   ierr = (*mat->ops->lufactor)(mat,row,col,f); CHKERRQ(ierr);
1135   PLogEventEnd(MAT_LUFactor,mat,row,col,0);
1136   PetscFunctionReturn(0);
1137 }
1138 
1139 #undef __FUNC__
1140 #define __FUNC__ "MatILUFactor"
1141 /*@
1142    MatILUFactor - Performs in-place ILU factorization of matrix.
1143 
1144    Collective on Mat
1145 
1146    Input Parameters:
1147 +  mat - the matrix
1148 .  row - row permutation
1149 .  col - column permutation
1150 -  info - structure containing
1151 $      levels - number of levels of fill.
1152 $      expected fill - as ratio of original fill.
1153 $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
1154                 missing diagonal entries)
1155 
1156    Notes:
1157    Probably really in-place only when level of fill is zero, otherwise allocates
1158    new space to store factored matrix and deletes previous memory.
1159 
1160 .keywords: matrix, factor, ILU, in-place
1161 
1162 .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
1163 @*/
1164 int MatILUFactor(Mat mat,IS row,IS col,MatILUInfo *info)
1165 {
1166   int ierr;
1167 
1168   PetscFunctionBegin;
1169   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1170   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square");
1171   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1172   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1173   if (!mat->ops->ilufactor) SETERRQ(PETSC_ERR_SUP,0,"");
1174 
1175   PLogEventBegin(MAT_ILUFactor,mat,row,col,0);
1176   ierr = (*mat->ops->ilufactor)(mat,row,col,info); CHKERRQ(ierr);
1177   PLogEventEnd(MAT_ILUFactor,mat,row,col,0);
1178   PetscFunctionReturn(0);
1179 }
1180 
1181 #undef __FUNC__
1182 #define __FUNC__ "MatLUFactorSymbolic"
1183 /*@
1184    MatLUFactorSymbolic - Performs symbolic LU factorization of matrix.
1185    Call this routine before calling MatLUFactorNumeric().
1186 
1187    Collective on Mat
1188 
1189    Input Parameters:
1190 +  mat - the matrix
1191 .  row, col - row and column permutations
1192 -  f - expected fill as ratio of the original number of nonzeros,
1193        for example 3.0; choosing this parameter well can result in
1194        more efficient use of time and space. Run with the option -log_info
1195        to determine an optimal value to use
1196 
1197    Output Parameter:
1198 .  fact - new matrix that has been symbolically factored
1199 
1200    Notes:
1201    See the file ${PETSC_DIR}/Performance for additional information about
1202    choosing the fill factor for better efficiency.
1203 
1204 .keywords: matrix, factor, LU, symbolic, fill
1205 
1206 .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor()
1207 @*/
1208 int MatLUFactorSymbolic(Mat mat,IS row,IS col,double f,Mat *fact)
1209 {
1210   int ierr;
1211 
1212   PetscFunctionBegin;
1213   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1214   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square");
1215   PetscValidPointer(fact);
1216   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1217   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1218   if (!mat->ops->lufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"");
1219 
1220   PLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);
1221   ierr = (*mat->ops->lufactorsymbolic)(mat,row,col,f,fact); CHKERRQ(ierr);
1222   PLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);
1223   PetscFunctionReturn(0);
1224 }
1225 
1226 #undef __FUNC__
1227 #define __FUNC__ "MatLUFactorNumeric"
1228 /*@
1229    MatLUFactorNumeric - Performs numeric LU factorization of a matrix.
1230    Call this routine after first calling MatLUFactorSymbolic().
1231 
1232    Collective on Mat
1233 
1234    Input Parameters:
1235 +  mat - the matrix
1236 -  row, col - row and  column permutations
1237 
1238    Output Parameters:
1239 .  fact - symbolically factored matrix that must have been generated
1240           by MatLUFactorSymbolic()
1241 
1242    Notes:
1243    See MatLUFactor() for in-place factorization.  See
1244    MatCholeskyFactorNumeric() for the symmetric, positive definite case.
1245 
1246 .keywords: matrix, factor, LU, numeric
1247 
1248 .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
1249 @*/
1250 int MatLUFactorNumeric(Mat mat,Mat *fact)
1251 {
1252   int ierr,flg;
1253 
1254   PetscFunctionBegin;
1255   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1256   PetscValidPointer(fact);
1257   PetscValidHeaderSpecific(*fact,MAT_COOKIE);
1258   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1259   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1260     SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dimensions are different %d should = %d %d should = %d",
1261             mat->M,(*fact)->M,mat->N,(*fact)->N);
1262   }
1263   if (!(*fact)->ops->lufactornumeric) SETERRQ(PETSC_ERR_SUP,0,"");
1264 
1265   PLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0);
1266   ierr = (*(*fact)->ops->lufactornumeric)(mat,fact); CHKERRQ(ierr);
1267   PLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0);
1268   ierr = OptionsHasName(PETSC_NULL,"-mat_view_draw",&flg); CHKERRQ(ierr);
1269   if (flg) {
1270     ierr = OptionsHasName(0,"-mat_view_contour",&flg); CHKERRQ(ierr);
1271     if (flg) {
1272       ViewerPushFormat(VIEWER_DRAW_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr);
1273     }
1274     ierr = MatView(*fact,VIEWER_DRAW_(mat->comm)); CHKERRQ(ierr);
1275     ierr = ViewerFlush(VIEWER_DRAW_(mat->comm)); CHKERRQ(ierr);
1276     if (flg) {
1277       ViewerPopFormat(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
1278     }
1279   }
1280   PetscFunctionReturn(0);
1281 }
1282 
1283 #undef __FUNC__
1284 #define __FUNC__ "MatCholeskyFactor"
1285 /*@
1286    MatCholeskyFactor - Performs in-place Cholesky factorization of a
1287    symmetric matrix.
1288 
1289    Collective on Mat
1290 
1291    Input Parameters:
1292 +  mat - the matrix
1293 .  perm - row and column permutations
1294 -  f - expected fill as ratio of original fill
1295 
1296    Notes:
1297    See MatLUFactor() for the nonsymmetric case.  See also
1298    MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().
1299 
1300 .keywords: matrix, factor, in-place, Cholesky
1301 
1302 .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
1303 @*/
1304 int MatCholeskyFactor(Mat mat,IS perm,double f)
1305 {
1306   int ierr;
1307 
1308   PetscFunctionBegin;
1309   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1310   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Matrix must be square");
1311   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1312   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1313   if (!mat->ops->choleskyfactor) SETERRQ(PETSC_ERR_SUP,0,"");
1314 
1315   PLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);
1316   ierr = (*mat->ops->choleskyfactor)(mat,perm,f); CHKERRQ(ierr);
1317   PLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);
1318   PetscFunctionReturn(0);
1319 }
1320 
1321 #undef __FUNC__
1322 #define __FUNC__ "MatCholeskyFactorSymbolic"
1323 /*@
1324    MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
1325    of a symmetric matrix.
1326 
1327    Collective on Mat
1328 
1329    Input Parameters:
1330 +  mat - the matrix
1331 .  perm - row and column permutations
1332 -  f - expected fill as ratio of original
1333 
1334    Output Parameter:
1335 .  fact - the factored matrix
1336 
1337    Notes:
1338    See MatLUFactorSymbolic() for the nonsymmetric case.  See also
1339    MatCholeskyFactor() and MatCholeskyFactorNumeric().
1340 
1341 .keywords: matrix, factor, factorization, symbolic, Cholesky
1342 
1343 .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
1344 @*/
1345 int MatCholeskyFactorSymbolic(Mat mat,IS perm,double f,Mat *fact)
1346 {
1347   int ierr;
1348 
1349   PetscFunctionBegin;
1350   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1351   PetscValidPointer(fact);
1352   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Matrix must be square");
1353   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1354   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1355   if (!mat->ops->choleskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"");
1356 
1357   PLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
1358   ierr = (*mat->ops->choleskyfactorsymbolic)(mat,perm,f,fact); CHKERRQ(ierr);
1359   PLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
1360   PetscFunctionReturn(0);
1361 }
1362 
1363 #undef __FUNC__
1364 #define __FUNC__ "MatCholeskyFactorNumeric"
1365 /*@
1366    MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
1367    of a symmetric matrix. Call this routine after first calling
1368    MatCholeskyFactorSymbolic().
1369 
1370    Collective on Mat
1371 
1372    Input Parameter:
1373 .  mat - the initial matrix
1374 
1375    Output Parameter:
1376 .  fact - the factored matrix
1377 
1378 .keywords: matrix, factor, numeric, Cholesky
1379 
1380 .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
1381 @*/
1382 int MatCholeskyFactorNumeric(Mat mat,Mat *fact)
1383 {
1384   int ierr;
1385 
1386   PetscFunctionBegin;
1387   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1388   PetscValidPointer(fact);
1389   if (!mat->ops->choleskyfactornumeric) SETERRQ(PETSC_ERR_SUP,0,"");
1390   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1391   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1392     SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dim %d should = %d %d should = %d",
1393             mat->M,(*fact)->M,mat->N,(*fact)->N);
1394   }
1395 
1396   PLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
1397   ierr = (*mat->ops->choleskyfactornumeric)(mat,fact); CHKERRQ(ierr);
1398   PLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
1399   PetscFunctionReturn(0);
1400 }
1401 
1402 /* ----------------------------------------------------------------*/
1403 #undef __FUNC__
1404 #define __FUNC__ "MatSolve"
1405 /*@
1406    MatSolve - Solves A x = b, given a factored matrix.
1407 
1408    Collective on Mat and Vec
1409 
1410    Input Parameters:
1411 +  mat - the factored matrix
1412 -  b - the right-hand-side vector
1413 
1414    Output Parameter:
1415 .  x - the result vector
1416 
1417    Notes:
1418    The vectors b and x cannot be the same.  I.e., one cannot
1419    call MatSolve(A,x,x).
1420 
1421 .keywords: matrix, linear system, solve, LU, Cholesky, triangular solve
1422 
1423 .seealso: MatSolveAdd(), MatSolveTrans(), MatSolveTransAdd()
1424 @*/
1425 int MatSolve(Mat mat,Vec b,Vec x)
1426 {
1427   int ierr;
1428 
1429   PetscFunctionBegin;
1430   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1431   PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE);
1432   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1433   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1434   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1435   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1436   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1437   if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0);
1438 
1439   if (!mat->ops->solve) SETERRQ(PETSC_ERR_SUP,0,"");
1440   PLogEventBegin(MAT_Solve,mat,b,x,0);
1441   ierr = (*mat->ops->solve)(mat,b,x); CHKERRQ(ierr);
1442   PLogEventEnd(MAT_Solve,mat,b,x,0);
1443   PetscFunctionReturn(0);
1444 }
1445 
1446 #undef __FUNC__
1447 #define __FUNC__ "MatForwardSolve"
1448 /* @
1449    MatForwardSolve - Solves L x = b, given a factored matrix, A = LU.
1450 
1451    Collective on Mat and Vec
1452 
1453    Input Parameters:
1454 +  mat - the factored matrix
1455 -  b - the right-hand-side vector
1456 
1457    Output Parameter:
1458 .  x - the result vector
1459 
1460    Notes:
1461    MatSolve() should be used for most applications, as it performs
1462    a forward solve followed by a backward solve.
1463 
1464    The vectors b and x cannot be the same.  I.e., one cannot
1465    call MatForwardSolve(A,x,x).
1466 
1467 .keywords: matrix, forward, LU, Cholesky, triangular solve
1468 
1469 .seealso: MatSolve(), MatBackwardSolve()
1470 @ */
1471 int MatForwardSolve(Mat mat,Vec b,Vec x)
1472 {
1473   int ierr;
1474 
1475   PetscFunctionBegin;
1476   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1477   PetscValidHeaderSpecific(b,VEC_COOKIE);  PetscValidHeaderSpecific(x,VEC_COOKIE);
1478   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1479   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1480   if (!mat->ops->forwardsolve) SETERRQ(PETSC_ERR_SUP,0,"");
1481   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1482   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1483   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1484 
1485   PLogEventBegin(MAT_ForwardSolve,mat,b,x,0);
1486   ierr = (*mat->ops->forwardsolve)(mat,b,x); CHKERRQ(ierr);
1487   PLogEventEnd(MAT_ForwardSolve,mat,b,x,0);
1488   PetscFunctionReturn(0);
1489 }
1490 
1491 #undef __FUNC__
1492 #define __FUNC__ "MatBackwardSolve"
1493 /* @
1494    MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.
1495 
1496    Collective on Mat and Vec
1497 
1498    Input Parameters:
1499 +  mat - the factored matrix
1500 -  b - the right-hand-side vector
1501 
1502    Output Parameter:
1503 .  x - the result vector
1504 
1505    Notes:
1506    MatSolve() should be used for most applications, as it performs
1507    a forward solve followed by a backward solve.
1508 
1509    The vectors b and x cannot be the same.  I.e., one cannot
1510    call MatBackwardSolve(A,x,x).
1511 
1512 .keywords: matrix, backward, LU, Cholesky, triangular solve
1513 
1514 .seealso: MatSolve(), MatForwardSolve()
1515 @ */
1516 int MatBackwardSolve(Mat mat,Vec b,Vec x)
1517 {
1518   int ierr;
1519 
1520   PetscFunctionBegin;
1521   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1522   PetscValidHeaderSpecific(b,VEC_COOKIE);  PetscValidHeaderSpecific(x,VEC_COOKIE);
1523   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1524   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1525   if (!mat->ops->backwardsolve) SETERRQ(PETSC_ERR_SUP,0,"");
1526   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1527   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1528   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1529 
1530   PLogEventBegin(MAT_BackwardSolve,mat,b,x,0);
1531   ierr = (*mat->ops->backwardsolve)(mat,b,x); CHKERRQ(ierr);
1532   PLogEventEnd(MAT_BackwardSolve,mat,b,x,0);
1533   PetscFunctionReturn(0);
1534 }
1535 
1536 #undef __FUNC__
1537 #define __FUNC__ "MatSolveAdd"
1538 /*@
1539    MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.
1540 
1541    Collective on Mat and Vec
1542 
1543    Input Parameters:
1544 +  mat - the factored matrix
1545 .  b - the right-hand-side vector
1546 -  y - the vector to be added to
1547 
1548    Output Parameter:
1549 .  x - the result vector
1550 
1551    Notes:
1552    The vectors b and x cannot be the same.  I.e., one cannot
1553    call MatSolveAdd(A,x,y,x).
1554 
1555 .keywords: matrix, linear system, solve, LU, Cholesky, add
1556 
1557 .seealso: MatSolve(), MatSolveTrans(), MatSolveTransAdd()
1558 @*/
1559 int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
1560 {
1561   Scalar one = 1.0;
1562   Vec    tmp;
1563   int    ierr;
1564 
1565   PetscFunctionBegin;
1566   PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE);
1567   PetscValidHeaderSpecific(b,VEC_COOKIE);  PetscValidHeaderSpecific(x,VEC_COOKIE);
1568   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1569   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1570   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1571   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1572   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
1573   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1574   if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim %d %d",x->n,y->n);
1575 
1576   PLogEventBegin(MAT_SolveAdd,mat,b,x,y);
1577   if (mat->ops->solveadd)  {
1578     ierr = (*mat->ops->solveadd)(mat,b,y,x); CHKERRQ(ierr);
1579   } else {
1580     /* do the solve then the add manually */
1581     if (x != y) {
1582       ierr = MatSolve(mat,b,x); CHKERRQ(ierr);
1583       ierr = VecAXPY(&one,y,x); CHKERRQ(ierr);
1584     } else {
1585       ierr = VecDuplicate(x,&tmp); CHKERRQ(ierr);
1586       PLogObjectParent(mat,tmp);
1587       ierr = VecCopy(x,tmp); CHKERRQ(ierr);
1588       ierr = MatSolve(mat,b,x); CHKERRQ(ierr);
1589       ierr = VecAXPY(&one,tmp,x); CHKERRQ(ierr);
1590       ierr = VecDestroy(tmp); CHKERRQ(ierr);
1591     }
1592   }
1593   PLogEventEnd(MAT_SolveAdd,mat,b,x,y);
1594   PetscFunctionReturn(0);
1595 }
1596 
1597 #undef __FUNC__
1598 #define __FUNC__ "MatSolveTrans"
1599 /*@
1600    MatSolveTrans - Solves A' x = b, given a factored matrix.
1601 
1602    Collective on Mat and Vec
1603 
1604    Input Parameters:
1605 +  mat - the factored matrix
1606 -  b - the right-hand-side vector
1607 
1608    Output Parameter:
1609 .  x - the result vector
1610 
1611    Notes:
1612    The vectors b and x cannot be the same.  I.e., one cannot
1613    call MatSolveTrans(A,x,x).
1614 
1615 .keywords: matrix, linear system, solve, LU, Cholesky, transpose
1616 
1617 .seealso: MatSolve(), MatSolveAdd(), MatSolveTransAdd()
1618 @*/
1619 int MatSolveTrans(Mat mat,Vec b,Vec x)
1620 {
1621   int ierr;
1622 
1623   PetscFunctionBegin;
1624   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1625   PetscValidHeaderSpecific(b,VEC_COOKIE);  PetscValidHeaderSpecific(x,VEC_COOKIE);
1626   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1627   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1628   if (!mat->ops->solvetrans) SETERRQ(PETSC_ERR_SUP,0,"");
1629   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
1630   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);
1631 
1632   PLogEventBegin(MAT_SolveTrans,mat,b,x,0);
1633   ierr = (*mat->ops->solvetrans)(mat,b,x); CHKERRQ(ierr);
1634   PLogEventEnd(MAT_SolveTrans,mat,b,x,0);
1635   PetscFunctionReturn(0);
1636 }
1637 
1638 #undef __FUNC__
1639 #define __FUNC__ "MatSolveTransAdd"
1640 /*@
1641    MatSolveTransAdd - Computes x = y + inv(trans(A)) b, given a
1642                       factored matrix.
1643 
1644    Collective on Mat and Vec
1645 
1646    Input Parameters:
1647 +  mat - the factored matrix
1648 .  b - the right-hand-side vector
1649 -  y - the vector to be added to
1650 
1651    Output Parameter:
1652 .  x - the result vector
1653 
1654    Notes:
1655    The vectors b and x cannot be the same.  I.e., one cannot
1656    call MatSolveTransAdd(A,x,y,x).
1657 
1658 .keywords: matrix, linear system, solve, LU, Cholesky, transpose, add
1659 
1660 .seealso: MatSolve(), MatSolveAdd(), MatSolveTrans()
1661 @*/
1662 int MatSolveTransAdd(Mat mat,Vec b,Vec y,Vec x)
1663 {
1664   Scalar one = 1.0;
1665   int    ierr;
1666   Vec    tmp;
1667 
1668   PetscFunctionBegin;
1669   PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE);
1670   PetscValidHeaderSpecific(b,VEC_COOKIE);  PetscValidHeaderSpecific(x,VEC_COOKIE);
1671   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors");
1672   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix");
1673   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
1674   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);
1675   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
1676   if (x->n != y->n)   SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim %d %d",x->n,y->n);
1677 
1678   PLogEventBegin(MAT_SolveTransAdd,mat,b,x,y);
1679   if (mat->ops->solvetransadd) {
1680     ierr = (*mat->ops->solvetransadd)(mat,b,y,x); CHKERRQ(ierr);
1681   } else {
1682     /* do the solve then the add manually */
1683     if (x != y) {
1684       ierr = MatSolveTrans(mat,b,x); CHKERRQ(ierr);
1685       ierr = VecAXPY(&one,y,x); CHKERRQ(ierr);
1686     } else {
1687       ierr = VecDuplicate(x,&tmp); CHKERRQ(ierr);
1688       PLogObjectParent(mat,tmp);
1689       ierr = VecCopy(x,tmp); CHKERRQ(ierr);
1690       ierr = MatSolveTrans(mat,b,x); CHKERRQ(ierr);
1691       ierr = VecAXPY(&one,tmp,x); CHKERRQ(ierr);
1692       ierr = VecDestroy(tmp); CHKERRQ(ierr);
1693     }
1694   }
1695   PLogEventEnd(MAT_SolveTransAdd,mat,b,x,y);
1696   PetscFunctionReturn(0);
1697 }
1698 /* ----------------------------------------------------------------*/
1699 
1700 #undef __FUNC__
1701 #define __FUNC__ "MatRelax"
1702 /*@
1703    MatRelax - Computes one relaxation sweep.
1704 
1705    Collective on Mat and Vec
1706 
1707    Input Parameters:
1708 +  mat - the matrix
1709 .  b - the right hand side
1710 .  omega - the relaxation factor
1711 .  flag - flag indicating the type of SOR (see below)
1712 .  shift -  diagonal shift
1713 -  its - the number of iterations
1714 
1715    Output Parameters:
1716 .  x - the solution (can contain an initial guess)
1717 
1718    SOR Flags:
1719 .     SOR_FORWARD_SWEEP - forward SOR
1720 .     SOR_BACKWARD_SWEEP - backward SOR
1721 .     SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
1722 .     SOR_LOCAL_FORWARD_SWEEP - local forward SOR
1723 .     SOR_LOCAL_BACKWARD_SWEEP - local forward SOR
1724 .     SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
1725 .     SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies
1726          upper/lower triangular part of matrix to
1727          vector (with omega)
1728 .     SOR_ZERO_INITIAL_GUESS - zero initial guess
1729 
1730    Notes:
1731    SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
1732    SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings
1733    on each processor.
1734 
1735    Application programmers will not generally use MatRelax() directly,
1736    but instead will employ the SLES/PC interface.
1737 
1738    Notes for Advanced Users:
1739    The flags are implemented as bitwise inclusive or operations.
1740    For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
1741    to specify a zero initial guess for SSOR.
1742 
1743 .keywords: matrix, relax, relaxation, sweep
1744 @*/
1745 int MatRelax(Mat mat,Vec b,double omega,MatSORType flag,double shift,
1746              int its,Vec x)
1747 {
1748   int ierr;
1749 
1750   PetscFunctionBegin;
1751   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1752   PetscValidHeaderSpecific(b,VEC_COOKIE);  PetscValidHeaderSpecific(x,VEC_COOKIE);
1753   if (!mat->ops->relax) SETERRQ(PETSC_ERR_SUP,0,"");
1754   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1755   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1756   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1757   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1758   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1759 
1760   PLogEventBegin(MAT_Relax,mat,b,x,0);
1761   ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,x); CHKERRQ(ierr);
1762   PLogEventEnd(MAT_Relax,mat,b,x,0);
1763   PetscFunctionReturn(0);
1764 }
1765 
1766 #undef __FUNC__
1767 #define __FUNC__ "MatCopy_Basic"
1768 /*
1769       Default matrix copy routine.
1770 */
1771 int MatCopy_Basic(Mat A,Mat B,MatStructure str)
1772 {
1773   int    ierr,i,rstart,rend,nz,*cwork;
1774   Scalar *vwork;
1775 
1776   PetscFunctionBegin;
1777   ierr = MatZeroEntries(B); CHKERRQ(ierr);
1778   ierr = MatGetOwnershipRange(A,&rstart,&rend); CHKERRQ(ierr);
1779   for (i=rstart; i<rend; i++) {
1780     ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr);
1781     ierr = MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr);
1782     ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr);
1783   }
1784   ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1785   ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1786   PetscFunctionReturn(0);
1787 }
1788 
1789 #undef __FUNC__
1790 #define __FUNC__ "MatCopy"
1791 /*@C
1792    MatCopy - Copys a matrix to another matrix.
1793 
1794    Collective on Mat
1795 
1796    Input Parameters:
1797 +  A - the matrix
1798 -  str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN
1799 
1800    Output Parameter:
1801 .  B - where the copy is put
1802 
1803    Notes:
1804    MatCopy() copies the matrix entries of a matrix to another existing
1805    matrix (after first zeroing the second matrix).  A related routine is
1806    MatConvert(), which first creates a new matrix and then copies the data.
1807 
1808 .keywords: matrix, copy, convert
1809 
1810 .seealso: MatConvert()
1811 @*/
1812 int MatCopy(Mat A,Mat B,MatStructure str)
1813 {
1814   int ierr;
1815 
1816   PetscFunctionBegin;
1817   PetscValidHeaderSpecific(A,MAT_COOKIE); PetscValidHeaderSpecific(B,MAT_COOKIE);
1818   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1819   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1820   if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat A,Mat B: global dim %d %d",A->M,B->M,
1821                                              A->N,B->N);
1822 
1823   PLogEventBegin(MAT_Copy,A,B,0,0);
1824   if (A->ops->copy) {
1825     ierr = (*A->ops->copy)(A,B,str); CHKERRQ(ierr);
1826   } else { /* generic conversion */
1827     ierr = MatCopy_Basic(A,B,str); CHKERRQ(ierr);
1828   }
1829   PLogEventEnd(MAT_Copy,A,B,0,0);
1830   PetscFunctionReturn(0);
1831 }
1832 
1833 static int MatConvertersSet = 0;
1834 static int (*MatConverters[MAX_MATRIX_TYPES][MAX_MATRIX_TYPES])(Mat,MatType,Mat*) =
1835            {{0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
1836             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
1837             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
1838             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
1839             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0},
1840             {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}};
1841 
1842 #undef __FUNC__
1843 #define __FUNC__ "MatConvertRegister"
1844 /*@C
1845     MatConvertRegister - Allows one to register a routine that converts between
1846         two matrix types.
1847 
1848     Not Collective
1849 
1850     Input Parameters:
1851 .   intype - the type of matrix (defined in include/mat.h), for example, MATSEQAIJ.
1852 .   outtype - new matrix type, or MATSAME
1853 
1854 .seealso: MatConvertRegisterAll()
1855 @*/
1856 int MatConvertRegister(MatType intype,MatType outtype,int (*converter)(Mat,MatType,Mat*))
1857 {
1858   PetscFunctionBegin;
1859   MatConverters[intype][outtype] = converter;
1860   MatConvertersSet               = 1;
1861   PetscFunctionReturn(0);
1862 }
1863 
1864 #undef __FUNC__
1865 #define __FUNC__ "MatConvert"
1866 /*@C
1867    MatConvert - Converts a matrix to another matrix, either of the same
1868    or different type.
1869 
1870    Collective on Mat
1871 
1872    Input Parameters:
1873 +  mat - the matrix
1874 -  newtype - new matrix type.  Use MATSAME to create a new matrix of the
1875    same type as the original matrix.
1876 
1877    Output Parameter:
1878 .  M - pointer to place new matrix
1879 
1880    Notes:
1881    MatConvert() first creates a new matrix and then copies the data from
1882    the first matrix.  A related routine is MatCopy(), which copies the matrix
1883    entries of one matrix to another already existing matrix context.
1884 
1885 .keywords: matrix, copy, convert
1886 
1887 .seealso: MatCopy(), MatDuplicate()
1888 @*/
1889 int MatConvert(Mat mat,MatType newtype,Mat *M)
1890 {
1891   int ierr;
1892 
1893   PetscFunctionBegin;
1894   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1895   PetscValidPointer(M);
1896   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1897   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1898 
1899   if (newtype > MAX_MATRIX_TYPES || newtype < -1) {
1900     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,1,"Not a valid matrix type");
1901   }
1902   *M  = 0;
1903 
1904   if (!MatConvertersSet) {
1905     ierr = MatLoadRegisterAll(); CHKERRQ(ierr);
1906   }
1907 
1908   PLogEventBegin(MAT_Convert,mat,0,0,0);
1909   if ((newtype == mat->type || newtype == MATSAME) && mat->ops->duplicate) {
1910     ierr = (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M); CHKERRQ(ierr);
1911   } else {
1912     if (!MatConvertersSet) {
1913       ierr = MatConvertRegisterAll(); CHKERRQ(ierr);
1914     }
1915     if (!MatConverters[mat->type][newtype]) {
1916       SETERRQ(PETSC_ERR_ARG_WRONG,1,"Invalid matrix type, or matrix converter not registered");
1917     }
1918     ierr = (*MatConverters[mat->type][newtype])(mat,newtype,M); CHKERRQ(ierr);
1919   }
1920   PLogEventEnd(MAT_Convert,mat,0,0,0);
1921   PetscFunctionReturn(0);
1922 }
1923 
1924 #undef __FUNC__
1925 #define __FUNC__ "MatDuplicate"
1926 /*@C
1927    MatDuplicate - Duplicates a matrix including the non-zero structure.
1928 
1929    Collective on Mat
1930 
1931    Input Parameters:
1932 +  mat - the matrix
1933 -  op - either MAT_DO_NO_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero
1934         values as well or not
1935 
1936    Output Parameter:
1937 .  M - pointer to place new matrix
1938 
1939 .keywords: matrix, copy, convert, duplicate
1940 
1941 .seealso: MatCopy(), MatConvert()
1942 @*/
1943 int MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
1944 {
1945   int ierr;
1946 
1947   PetscFunctionBegin;
1948   PetscValidHeaderSpecific(mat,MAT_COOKIE);
1949   PetscValidPointer(M);
1950   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1951   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
1952 
1953   *M  = 0;
1954   PLogEventBegin(MAT_Convert,mat,0,0,0);
1955   if (!mat->ops->duplicate) {
1956     SETERRQ(PETSC_ERR_SUP,1,"Not written for this matrix type");
1957   }
1958   ierr = (*mat->ops->duplicate)(mat,op,M); CHKERRQ(ierr);
1959   PLogEventEnd(MAT_Convert,mat,0,0,0);
1960   PetscFunctionReturn(0);
1961 }
1962 
1963 #undef __FUNC__
1964 #define __FUNC__ "MatGetDiagonal"
1965 /*@
1966    MatGetDiagonal - Gets the diagonal of a matrix.
1967 
1968    Collective on Mat and Vec
1969 
1970    Input Parameters:
1971 +  mat - the matrix
1972 -  v - the vector for storing the diagonal
1973 
1974    Output Parameter:
1975 .  v - the diagonal of the matrix
1976 
1977    Notes:
1978    For the SeqAIJ matrix format, this routine may also be called
1979    on a LU factored matrix; in that case it routines the reciprocal of
1980    the diagonal entries in U. It returns the entries permuted by the
1981    row and column permutation used during the symbolic factorization.
1982 
1983 .keywords: matrix, get, diagonal
1984 @*/
1985 int MatGetDiagonal(Mat mat,Vec v)
1986 {
1987   int ierr;
1988 
1989   PetscFunctionBegin;
1990   PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(v,VEC_COOKIE);
1991   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
1992   if (!mat->ops->getdiagonal) SETERRQ(PETSC_ERR_SUP,0,"");
1993   ierr = (*mat->ops->getdiagonal)(mat,v);CHKERRQ(ierr);
1994   PetscFunctionReturn(0);
1995 }
1996 
1997 #undef __FUNC__
1998 #define __FUNC__ "MatTranspose"
1999 /*@C
2000    MatTranspose - Computes an in-place or out-of-place transpose of a matrix.
2001 
2002    Collective on Mat
2003 
2004    Input Parameter:
2005 .  mat - the matrix to transpose
2006 
2007    Output Parameters:
2008 .  B - the transpose (or pass in PETSC_NULL for an in-place transpose)
2009 
2010 .keywords: matrix, transpose
2011 
2012 .seealso: MatMultTrans(), MatMultTransAdd()
2013 @*/
2014 int MatTranspose(Mat mat,Mat *B)
2015 {
2016   int ierr;
2017 
2018   PetscFunctionBegin;
2019   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2020   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2021   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2022   if (!mat->ops->transpose) SETERRQ(PETSC_ERR_SUP,0,"");
2023   ierr = (*mat->ops->transpose)(mat,B);CHKERRQ(ierr);
2024   PetscFunctionReturn(0);
2025 }
2026 
2027 #undef __FUNC__
2028 #define __FUNC__ "MatPermute"
2029 /*@C
2030    MatPermute - Creates a new matrix with rows and columns permuted from the
2031    original.
2032 
2033    Collective on Mat
2034 
2035    Input Parameters:
2036 +  mat - the matrix to permute
2037 .  row - row permutation
2038 -  col - column permutation
2039 
2040    Output Parameters:
2041 .  B - the permuted matrix
2042 
2043 .keywords: matrix, transpose
2044 
2045 .seealso: MatGetReordering()
2046 @*/
2047 int MatPermute(Mat mat,IS row,IS col,Mat *B)
2048 {
2049   int ierr;
2050 
2051   PetscFunctionBegin;
2052   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2053   PetscValidHeaderSpecific(row,IS_COOKIE);
2054   PetscValidHeaderSpecific(col,IS_COOKIE);
2055   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2056   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2057   if (!mat->ops->permute) SETERRQ(PETSC_ERR_SUP,0,"");
2058   ierr = (*mat->ops->permute)(mat,row,col,B);CHKERRQ(ierr);
2059   PetscFunctionReturn(0);
2060 }
2061 
2062 #undef __FUNC__
2063 #define __FUNC__ "MatEqual"
2064 /*@
2065    MatEqual - Compares two matrices.
2066 
2067    Collective on Mat
2068 
2069    Input Parameters:
2070 +  A - the first matrix
2071 -  B - the second matrix
2072 
2073    Output Parameter:
2074 .  flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.
2075 
2076 .keywords: matrix, equal, equivalent
2077 @*/
2078 int MatEqual(Mat A,Mat B,PetscTruth *flg)
2079 {
2080   int ierr;
2081 
2082   PetscFunctionBegin;
2083   PetscValidHeaderSpecific(A,MAT_COOKIE); PetscValidHeaderSpecific(B,MAT_COOKIE);
2084   PetscValidIntPointer(flg);
2085   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2086   if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2087   if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat A,Mat B: global dim %d %d %d %d",
2088                                              A->M,B->M,A->N,B->N);
2089   if (!A->ops->equal) SETERRQ(PETSC_ERR_SUP,0,"");
2090   ierr = (*A->ops->equal)(A,B,flg);CHKERRQ(ierr);
2091   PetscFunctionReturn(0);
2092 }
2093 
2094 #undef __FUNC__
2095 #define __FUNC__ "MatDiagonalScale"
2096 /*@
2097    MatDiagonalScale - Scales a matrix on the left and right by diagonal
2098    matrices that are stored as vectors.  Either of the two scaling
2099    matrices can be PETSC_NULL.
2100 
2101    Collective on Mat
2102 
2103    Input Parameters:
2104 +  mat - the matrix to be scaled
2105 .  l - the left scaling vector (or PETSC_NULL)
2106 -  r - the right scaling vector (or PETSC_NULL)
2107 
2108    Notes:
2109    MatDiagonalScale() computes A = LAR, where
2110    L = a diagonal matrix, R = a diagonal matrix
2111 
2112 .keywords: matrix, diagonal, scale
2113 
2114 .seealso: MatDiagonalScale()
2115 @*/
2116 int MatDiagonalScale(Mat mat,Vec l,Vec r)
2117 {
2118   int ierr;
2119 
2120   PetscFunctionBegin;
2121   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2122   if (!mat->ops->diagonalscale) SETERRQ(PETSC_ERR_SUP,0,"");
2123   if (l) PetscValidHeaderSpecific(l,VEC_COOKIE);
2124   if (r) PetscValidHeaderSpecific(r,VEC_COOKIE);
2125   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2126   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2127 
2128   PLogEventBegin(MAT_Scale,mat,0,0,0);
2129   ierr = (*mat->ops->diagonalscale)(mat,l,r); CHKERRQ(ierr);
2130   PLogEventEnd(MAT_Scale,mat,0,0,0);
2131   PetscFunctionReturn(0);
2132 }
2133 
2134 #undef __FUNC__
2135 #define __FUNC__ "MatScale"
2136 /*@
2137     MatScale - Scales all elements of a matrix by a given number.
2138 
2139     Collective on Mat
2140 
2141     Input Parameters:
2142 +   mat - the matrix to be scaled
2143 -   a  - the scaling value
2144 
2145     Output Parameter:
2146 .   mat - the scaled matrix
2147 
2148 .keywords: matrix, scale
2149 
2150 .seealso: MatDiagonalScale()
2151 @*/
2152 int MatScale(Scalar *a,Mat mat)
2153 {
2154   int ierr;
2155 
2156   PetscFunctionBegin;
2157   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2158   PetscValidScalarPointer(a);
2159   if (!mat->ops->scale) SETERRQ(PETSC_ERR_SUP,0,"");
2160   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2161   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2162 
2163   PLogEventBegin(MAT_Scale,mat,0,0,0);
2164   ierr = (*mat->ops->scale)(a,mat); CHKERRQ(ierr);
2165   PLogEventEnd(MAT_Scale,mat,0,0,0);
2166   PetscFunctionReturn(0);
2167 }
2168 
2169 #undef __FUNC__
2170 #define __FUNC__ "MatNorm"
2171 /*@
2172    MatNorm - Calculates various norms of a matrix.
2173 
2174    Collective on Mat
2175 
2176    Input Parameters:
2177 +  mat - the matrix
2178 -  type - the type of norm, NORM_1, NORM_2, NORM_FROBENIUS, NORM_INFINITY
2179 
2180    Output Parameters:
2181 .  norm - the resulting norm
2182 
2183 .keywords: matrix, norm, Frobenius
2184 @*/
2185 int MatNorm(Mat mat,NormType type,double *norm)
2186 {
2187   int ierr;
2188 
2189   PetscFunctionBegin;
2190   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2191   PetscValidScalarPointer(norm);
2192 
2193   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2194   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2195   if (!mat->ops->norm) SETERRQ(PETSC_ERR_SUP,0,"Not for this matrix type");
2196   ierr = (*mat->ops->norm)(mat,type,norm);CHKERRQ(ierr);
2197   PetscFunctionReturn(0);
2198 }
2199 
2200 /*
2201      This variable is used to prevent counting of MatAssemblyBegin() that
2202    are called from within a MatAssemblyEnd().
2203 */
2204 static int MatAssemblyEnd_InUse = 0;
2205 #undef __FUNC__
2206 #define __FUNC__ "MatAssemblyBegin"
2207 /*@
2208    MatAssemblyBegin - Begins assembling the matrix.  This routine should
2209    be called after completing all calls to MatSetValues().
2210 
2211    Collective on Mat
2212 
2213    Input Parameters:
2214 +  mat - the matrix
2215 -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
2216 
2217    Notes:
2218    MatSetValues() generally caches the values.  The matrix is ready to
2219    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
2220    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
2221    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
2222    using the matrix.
2223 
2224    Level: beginner
2225 
2226 .keywords: matrix, assembly, assemble, begin
2227 
2228 .seealso: MatAssemblyEnd(), MatSetValues()
2229 @*/
2230 int MatAssemblyBegin(Mat mat,MatAssemblyType type)
2231 {
2232   int ierr;
2233 
2234   PetscFunctionBegin;
2235   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2236   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix.\n did you forget to call MatSetUnfactored()?");
2237   if (mat->assembled) {
2238     mat->was_assembled = PETSC_TRUE;
2239     mat->assembled     = PETSC_FALSE;
2240   }
2241   if (!MatAssemblyEnd_InUse) {
2242     PLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);
2243     if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);}
2244     PLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);
2245   } else {
2246     if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);}
2247   }
2248   PetscFunctionReturn(0);
2249 }
2250 
2251 
2252 #undef __FUNC__
2253 #define __FUNC__ "MatView_Private"
2254 /*
2255     Processes command line options to determine if/how a matrix
2256   is to be viewed. Called by MatAssemblyEnd() and MatLoad().
2257 */
2258 int MatView_Private(Mat mat)
2259 {
2260   int ierr,flg;
2261 
2262   PetscFunctionBegin;
2263   ierr = OptionsHasName(PETSC_NULL,"-mat_view_info",&flg); CHKERRQ(ierr);
2264   if (flg) {
2265     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO,0);CHKERRQ(ierr);
2266     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr);
2267     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2268   }
2269   ierr = OptionsHasName(PETSC_NULL,"-mat_view_info_detailed",&flg);CHKERRQ(ierr);
2270   if (flg) {
2271     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO_LONG,0);CHKERRQ(ierr);
2272     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr);
2273     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2274   }
2275   ierr = OptionsHasName(PETSC_NULL,"-mat_view",&flg); CHKERRQ(ierr);
2276   if (flg) {
2277     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr);
2278   }
2279   ierr = OptionsHasName(PETSC_NULL,"-mat_view_matlab",&flg); CHKERRQ(ierr);
2280   if (flg) {
2281     ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_MATLAB,"M");CHKERRQ(ierr);
2282     ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr);
2283     ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr);
2284   }
2285   ierr = OptionsHasName(PETSC_NULL,"-mat_view_draw",&flg); CHKERRQ(ierr);
2286   if (flg) {
2287     ierr = OptionsHasName(0,"-mat_view_contour",&flg); CHKERRQ(ierr);
2288     if (flg) {
2289       ViewerPushFormat(VIEWER_DRAW_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr);
2290     }
2291     ierr = MatView(mat,VIEWER_DRAW_(mat->comm)); CHKERRQ(ierr);
2292     ierr = ViewerFlush(VIEWER_DRAW_(mat->comm)); CHKERRQ(ierr);
2293     if (flg) {
2294       ViewerPopFormat(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr);
2295     }
2296   }
2297   ierr = OptionsHasName(PETSC_NULL,"-mat_view_socket",&flg); CHKERRQ(ierr);
2298   if (flg) {
2299     ierr = MatView(mat,VIEWER_SOCKET_(mat->comm)); CHKERRQ(ierr);
2300     ierr = ViewerFlush(VIEWER_SOCKET_(mat->comm)); CHKERRQ(ierr);
2301   }
2302   PetscFunctionReturn(0);
2303 }
2304 
2305 #undef __FUNC__
2306 #define __FUNC__ "MatAssemblyEnd"
2307 /*@
2308    MatAssemblyEnd - Completes assembling the matrix.  This routine should
2309    be called after MatAssemblyBegin().
2310 
2311    Collective on Mat
2312 
2313    Input Parameters:
2314 +  mat - the matrix
2315 -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
2316 
2317    Options Database Keys:
2318 +  -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
2319 .  -mat_view_info_detailed - Prints more detailed info
2320 .  -mat_view - Prints matrix in ASCII format
2321 .  -mat_view_matlab - Prints matrix in Matlab format
2322 .  -mat_view_draw - Draws nonzero structure of matrix, using MatView() and DrawOpenX().
2323 .  -display <name> - Sets display name (default is host)
2324 -  -draw_pause <sec> - Sets number of seconds to pause after display
2325 
2326    Notes:
2327    MatSetValues() generally caches the values.  The matrix is ready to
2328    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
2329    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
2330    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
2331    using the matrix.
2332 
2333    Level: beginner
2334 
2335 .keywords: matrix, assembly, assemble, end
2336 
2337 .seealso: MatAssemblyBegin(), MatSetValues(), DrawOpenX(), MatView()
2338 @*/
2339 int MatAssemblyEnd(Mat mat,MatAssemblyType type)
2340 {
2341   int        ierr;
2342   static int inassm = 0;
2343 
2344   PetscFunctionBegin;
2345   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2346 
2347   inassm++;
2348   MatAssemblyEnd_InUse++;
2349   if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
2350     PLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);
2351     if (mat->ops->assemblyend) {
2352       ierr = (*mat->ops->assemblyend)(mat,type); CHKERRQ(ierr);
2353     }
2354     PLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);
2355   } else {
2356     if (mat->ops->assemblyend) {
2357       ierr = (*mat->ops->assemblyend)(mat,type); CHKERRQ(ierr);
2358     }
2359   }
2360 
2361   /* Flush assembly is not a true assembly */
2362   if (type != MAT_FLUSH_ASSEMBLY) {
2363     mat->assembled  = PETSC_TRUE; mat->num_ass++;
2364   }
2365   mat->insertmode = NOT_SET_VALUES;
2366   MatAssemblyEnd_InUse--;
2367 
2368   if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) {
2369     ierr = MatView_Private(mat); CHKERRQ(ierr);
2370   }
2371   inassm--;
2372   PetscFunctionReturn(0);
2373 }
2374 
2375 
2376 #undef __FUNC__
2377 #define __FUNC__ "MatCompress"
2378 /*@
2379    MatCompress - Tries to store the matrix in as little space as
2380    possible.  May fail if memory is already fully used, since it
2381    tries to allocate new space.
2382 
2383    Collective on Mat
2384 
2385    Input Parameters:
2386 .  mat - the matrix
2387 
2388 .keywords: matrix, compress
2389 @*/
2390 int MatCompress(Mat mat)
2391 {
2392   int ierr;
2393 
2394   PetscFunctionBegin;
2395   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2396   if (mat->ops->compress) {ierr = (*mat->ops->compress)(mat);CHKERRQ(ierr);}
2397   PetscFunctionReturn(0);
2398 }
2399 
2400 #undef __FUNC__
2401 #define __FUNC__ "MatSetOption"
2402 /*@
2403    MatSetOption - Sets a parameter option for a matrix. Some options
2404    may be specific to certain storage formats.  Some options
2405    determine how values will be inserted (or added). Sorted,
2406    row-oriented input will generally assemble the fastest. The default
2407    is row-oriented, nonsorted input.
2408 
2409    Collective on Mat
2410 
2411    Input Parameters:
2412 +  mat - the matrix
2413 -  option - the option, one of those listed below (and possibly others),
2414              e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR
2415 
2416    Options Describing Matrix Structure:
2417 +    MAT_SYMMETRIC - symmetric in terms of both structure and value
2418 -    MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure
2419 
2420    Options For Use with MatSetValues():
2421    Insert a logically dense subblock, which can be
2422 +    MAT_ROW_ORIENTED - row-oriented
2423 .    MAT_COLUMN_ORIENTED - column-oriented
2424 .    MAT_ROWS_SORTED - sorted by row
2425 .    MAT_ROWS_UNSORTED - not sorted by row
2426 .    MAT_COLUMNS_SORTED - sorted by column
2427 -    MAT_COLUMNS_UNSORTED - not sorted by column
2428 
2429    Not these options reflect the data you pass in with MatSetValues(); it has
2430    nothing to do with how the data is stored internally in the matrix
2431    data structure.
2432 
2433    When (re)assembling a matrix, we can restrict the input for
2434    efficiency/debugging purposes.  These options include
2435 +    MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be
2436         allowed if they generate a new nonzero
2437 .    MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed
2438 .    MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if
2439          they generate a nonzero in a new diagonal (for block diagonal format only)
2440 .    MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only)
2441 .    MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries
2442 .    MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry
2443 -    MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly
2444 
2445    Notes:
2446    Some options are relevant only for particular matrix types and
2447    are thus ignored by others.  Other options are not supported by
2448    certain matrix types and will generate an error message if set.
2449 
2450    If using a Fortran 77 module to compute a matrix, one may need to
2451    use the column-oriented option (or convert to the row-oriented
2452    format).
2453 
2454    MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion
2455    that would generate a new entry in the nonzero structure is instead
2456    ignored.  Thus, if memory has not alredy been allocated for this particular
2457    data, then the insertion is ignored. For dense matrices, in which
2458    the entire array is allocated, no entries are ever ignored.
2459 
2460    MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion
2461    that would generate a new entry in the nonzero structure instead produces
2462    an error. (Currently supported for AIJ and BAIJ formats only.)
2463    This is a useful flag when using SAME_NONZERO_PATTERN in calling
2464    SLESSetOperators() to ensure that the nonzero pattern truely does
2465    remain unchanged.
2466 
2467    MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion
2468    that would generate a new entry that has not been preallocated will
2469    instead produce an error. (Currently supported for AIJ and BAIJ formats
2470    only.) This is a useful flag when debugging matrix memory preallocation.
2471 
2472    MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for
2473    other processors should be dropped, rather than stashed.
2474    This is useful if you know that the "owning" processor is also
2475    always generating the correct matrix entries, so that PETSc need
2476    not transfer duplicate entries generated on another processor.
2477 
2478    MAT_USE_HASH_TABLE indicates that a hash table be used to improve the
2479    searches during matrix assembly. When this flag is set, the hash table
2480    is created during the first Matrix Assembly. This hash table is
2481    used the next time through, during MatSetVaules()/MatSetVaulesBlocked()
2482    to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag
2483    should be used with MAT_USE_HASH_TABLE flag. This option is currently
2484    supported by MATMPIBAIJ format only.
2485 
2486 .keywords: matrix, option, row-oriented, column-oriented, sorted, nonzero
2487 @*/
2488 int MatSetOption(Mat mat,MatOption op)
2489 {
2490   int ierr;
2491 
2492   PetscFunctionBegin;
2493   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2494   if (mat->ops->setoption) {ierr = (*mat->ops->setoption)(mat,op);CHKERRQ(ierr);}
2495   PetscFunctionReturn(0);
2496 }
2497 
2498 #undef __FUNC__
2499 #define __FUNC__ "MatZeroEntries"
2500 /*@
2501    MatZeroEntries - Zeros all entries of a matrix.  For sparse matrices
2502    this routine retains the old nonzero structure.
2503 
2504    Collective on Mat
2505 
2506    Input Parameters:
2507 .  mat - the matrix
2508 
2509 .keywords: matrix, zero, entries
2510 
2511 .seealso: MatZeroRows()
2512 @*/
2513 int MatZeroEntries(Mat mat)
2514 {
2515   int ierr;
2516 
2517   PetscFunctionBegin;
2518   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2519   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2520   if (!mat->ops->zeroentries) SETERRQ(PETSC_ERR_SUP,0,"");
2521 
2522   PLogEventBegin(MAT_ZeroEntries,mat,0,0,0);
2523   ierr = (*mat->ops->zeroentries)(mat); CHKERRQ(ierr);
2524   PLogEventEnd(MAT_ZeroEntries,mat,0,0,0);
2525   PetscFunctionReturn(0);
2526 }
2527 
2528 #undef __FUNC__
2529 #define __FUNC__ "MatZeroRows"
2530 /*@
2531    MatZeroRows - Zeros all entries (except possibly the main diagonal)
2532    of a set of rows of a matrix.
2533 
2534    Collective on Mat
2535 
2536    Input Parameters:
2537 +  mat - the matrix
2538 .  is - index set of rows to remove
2539 -  diag - pointer to value put in all diagonals of eliminated rows.
2540           Note that diag is not a pointer to an array, but merely a
2541           pointer to a single value.
2542 
2543    Notes:
2544    For the AIJ matrix formats this removes the old nonzero structure,
2545    but does not release memory.  For the dense and block diagonal
2546    formats this does not alter the nonzero structure.
2547 
2548    The user can set a value in the diagonal entry (or for the AIJ and
2549    row formats can optionally remove the main diagonal entry from the
2550    nonzero structure as well, by passing a null pointer as the final
2551    argument).
2552 
2553    For the parallel case, all processes that share the matrix (i.e.,
2554    those in the communicator used for matrix creation) MUST call this
2555    routine, regardless of whether any rows being zeroed are owned by
2556    them.
2557 
2558 .keywords: matrix, zero, rows, boundary conditions
2559 
2560 .seealso: MatZeroEntries(),
2561 @*/
2562 int MatZeroRows(Mat mat,IS is, Scalar *diag)
2563 {
2564   int ierr;
2565 
2566   PetscFunctionBegin;
2567   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2568   PetscValidHeaderSpecific(is,IS_COOKIE);
2569   if (diag) PetscValidScalarPointer(diag);
2570   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2571   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2572   if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,0,"");
2573 
2574   ierr = (*mat->ops->zerorows)(mat,is,diag); CHKERRQ(ierr);
2575   ierr = MatView_Private(mat); CHKERRQ(ierr);
2576   PetscFunctionReturn(0);
2577 }
2578 
2579 #undef __FUNC__
2580 #define __FUNC__ "MatZeroRowsLocal"
2581 /*@
2582    MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
2583    of a set of rows of a matrix; using local numbering of rows.
2584 
2585    Collective on Mat
2586 
2587    Input Parameters:
2588 +  mat - the matrix
2589 .  is - index set of rows to remove
2590 -  diag - pointer to value put in all diagonals of eliminated rows.
2591           Note that diag is not a pointer to an array, but merely a
2592           pointer to a single value.
2593 
2594    Notes:
2595    For the AIJ matrix formats this removes the old nonzero structure,
2596    but does not release memory.  For the dense and block diagonal
2597    formats this does not alter the nonzero structure.
2598 
2599    The user can set a value in the diagonal entry (or for the AIJ and
2600    row formats can optionally remove the main diagonal entry from the
2601    nonzero structure as well, by passing a null pointer as the final
2602    argument).
2603 
2604 .keywords: matrix, zero, rows, boundary conditions
2605 
2606 .seealso: MatZeroEntries(),
2607 @*/
2608 int MatZeroRowsLocal(Mat mat,IS is, Scalar *diag)
2609 {
2610   int ierr;
2611   IS  newis;
2612 
2613   PetscFunctionBegin;
2614   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2615   PetscValidHeaderSpecific(is,IS_COOKIE);
2616   if (diag) PetscValidScalarPointer(diag);
2617   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2618   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2619   if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,0,"");
2620 
2621   ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis); CHKERRQ(ierr);
2622   ierr =  (*mat->ops->zerorows)(mat,newis,diag); CHKERRQ(ierr);
2623   ierr = ISDestroy(newis);
2624   PetscFunctionReturn(0);
2625 }
2626 
2627 #undef __FUNC__
2628 #define __FUNC__ "MatGetSize"
2629 /*@
2630    MatGetSize - Returns the numbers of rows and columns in a matrix.
2631 
2632    Not Collective
2633 
2634    Input Parameter:
2635 .  mat - the matrix
2636 
2637    Output Parameters:
2638 +  m - the number of global rows
2639 -  n - the number of global columns
2640 
2641 .keywords: matrix, dimension, size, rows, columns, global, get
2642 
2643 .seealso: MatGetLocalSize()
2644 @*/
2645 int MatGetSize(Mat mat,int *m,int* n)
2646 {
2647   int ierr;
2648 
2649   PetscFunctionBegin;
2650   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2651   ierr = (*mat->ops->getsize)(mat,m,n);CHKERRQ(ierr);
2652   PetscFunctionReturn(0);
2653 }
2654 
2655 #undef __FUNC__
2656 #define __FUNC__ "MatGetLocalSize"
2657 /*@
2658    MatGetLocalSize - Returns the number of rows and columns in a matrix
2659    stored locally.  This information may be implementation dependent, so
2660    use with care.
2661 
2662    Not Collective
2663 
2664    Input Parameters:
2665 .  mat - the matrix
2666 
2667    Output Parameters:
2668 +  m - the number of local rows
2669 -  n - the number of local columns
2670 
2671 .keywords: matrix, dimension, size, local, rows, columns, get
2672 
2673 .seealso: MatGetSize()
2674 @*/
2675 int MatGetLocalSize(Mat mat,int *m,int* n)
2676 {
2677   int ierr;
2678 
2679   PetscFunctionBegin;
2680   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2681   ierr = (*mat->ops->getlocalsize)(mat,m,n);CHKERRQ(ierr);
2682   PetscFunctionReturn(0);
2683 }
2684 
2685 #undef __FUNC__
2686 #define __FUNC__ "MatGetOwnershipRange"
2687 /*@
2688    MatGetOwnershipRange - Returns the range of matrix rows owned by
2689    this processor, assuming that the matrix is laid out with the first
2690    n1 rows on the first processor, the next n2 rows on the second, etc.
2691    For certain parallel layouts this range may not be well defined.
2692 
2693    Not Collective
2694 
2695    Input Parameters:
2696 .  mat - the matrix
2697 
2698    Output Parameters:
2699 +  m - the global index of the first local row
2700 -  n - one more than the global index of the last local row
2701 
2702 .keywords: matrix, get, range, ownership
2703 @*/
2704 int MatGetOwnershipRange(Mat mat,int *m,int* n)
2705 {
2706   int ierr;
2707 
2708   PetscFunctionBegin;
2709   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2710   PetscValidIntPointer(m);
2711   PetscValidIntPointer(n);
2712   if (!mat->ops->getownershiprange) SETERRQ(PETSC_ERR_SUP,0,"");
2713   ierr = (*mat->ops->getownershiprange)(mat,m,n);CHKERRQ(ierr);
2714   PetscFunctionReturn(0);
2715 }
2716 
2717 #undef __FUNC__
2718 #define __FUNC__ "MatILUFactorSymbolic"
2719 /*@
2720    MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
2721    Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric()
2722    to complete the factorization.
2723 
2724    Collective on Mat
2725 
2726    Input Parameters:
2727 +  mat - the matrix
2728 .  row - row permutation
2729 .  column - column permutation
2730 -  info - structure containing
2731 $      levels - number of levels of fill.
2732 $      expected fill - as ratio of original fill.
2733 $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
2734                 missing diagonal entries)
2735 
2736    Output Parameters:
2737 .  fact - new matrix that has been symbolically factored
2738 
2739    Notes:
2740    See the file ${PETSC_DIR}/Performace for additional information about
2741    choosing the fill factor for better efficiency.
2742 
2743 .keywords: matrix, factor, incomplete, ILU, symbolic, fill
2744 
2745 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric()
2746 @*/
2747 int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatILUInfo *info,Mat *fact)
2748 {
2749   int ierr;
2750 
2751   PetscFunctionBegin;
2752   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2753   PetscValidPointer(fact);
2754   if (info && info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Levels of fill negative %d",info->levels);
2755   if (info && info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",info->fill);
2756   if (!mat->ops->ilufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Only MatCreateMPIRowbs() matrices support parallel ILU");
2757   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2758   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2759 
2760   PLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);
2761   ierr = (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact); CHKERRQ(ierr);
2762   PLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);
2763   PetscFunctionReturn(0);
2764 }
2765 
2766 #undef __FUNC__
2767 #define __FUNC__ "MatIncompleteCholeskyFactorSymbolic"
2768 /*@
2769    MatIncompleteCholeskyFactorSymbolic - Performs symbolic incomplete
2770    Cholesky factorization for a symmetric matrix.  Use
2771    MatCholeskyFactorNumeric() to complete the factorization.
2772 
2773    Collective on Mat
2774 
2775    Input Parameters:
2776 +  mat - the matrix
2777 .  perm - row and column permutation
2778 .  fill - levels of fill
2779 -  f - expected fill as ratio of original fill
2780 
2781    Output Parameter:
2782 .  fact - the factored matrix
2783 
2784    Note:  Currently only no-fill factorization is supported.
2785 
2786 .keywords: matrix, factor, incomplete, ICC, Cholesky, symbolic, fill
2787 
2788 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor()
2789 @*/
2790 int MatIncompleteCholeskyFactorSymbolic(Mat mat,IS perm,double f,int fill,Mat *fact)
2791 {
2792   int ierr;
2793 
2794   PetscFunctionBegin;
2795   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2796   PetscValidPointer(fact);
2797   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
2798   if (fill < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Fill negative %d",fill);
2799   if (f < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",f);
2800   if (!mat->ops->incompletecholeskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support ICC in parallel");
2801   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2802 
2803   PLogEventBegin(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);
2804   ierr = (*mat->ops->incompletecholeskyfactorsymbolic)(mat,perm,f,fill,fact);CHKERRQ(ierr);
2805   PLogEventEnd(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0);
2806   PetscFunctionReturn(0);
2807 }
2808 
2809 #undef __FUNC__
2810 #define __FUNC__ "MatGetArray"
2811 /*@C
2812    MatGetArray - Returns a pointer to the element values in the matrix.
2813    The result of this routine is dependent on the underlying matrix data
2814    structure, and may not even work for certain matrix types.  You MUST
2815    call MatRestoreArray() when you no longer need to access the array.
2816 
2817    Not Collective
2818 
2819    Input Parameter:
2820 .  mat - the matrix
2821 
2822    Output Parameter:
2823 .  v - the location of the values
2824 
2825    Currently returns an array only for the dense formats, giving access to
2826    the local portion of the matrix in the usual Fortran column-oriented format.
2827 
2828    Fortran Note:
2829    This routine is used differently from Fortran
2830 $    Mat         mat
2831 $    Scalar      mat_array(1)
2832 $    PetscOffset i_mat
2833 $    int         ierr
2834 $       call MatGetArray(mat,mat_array,i_mat,ierr)
2835 $
2836 $   Access first local entry in matrix with
2837 $      value = mat_array(i_mat + 1)
2838 $   (note here array is treated as one dimensional)
2839 $      ...... other code
2840 $       call MatRestoreArray(mat,mat_array,i_mat,ierr)
2841 
2842    See the Fortran chapter of the users manual and
2843    petsc/src/mat/examples/tests for details.
2844 
2845 .keywords: matrix, array, elements, values
2846 
2847 .seealso: MatRestoreArray()
2848 @*/
2849 int MatGetArray(Mat mat,Scalar **v)
2850 {
2851   int ierr;
2852 
2853   PetscFunctionBegin;
2854   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2855   PetscValidPointer(v);
2856   if (!mat->ops->getarray) SETERRQ(PETSC_ERR_SUP,0,"");
2857   ierr = (*mat->ops->getarray)(mat,v);CHKERRQ(ierr);
2858   PetscFunctionReturn(0);
2859 }
2860 
2861 #undef __FUNC__
2862 #define __FUNC__ "MatRestoreArray"
2863 /*@C
2864    MatRestoreArray - Restores the matrix after MatGetArray() has been called.
2865 
2866    Not Collective
2867 
2868    Input Parameter:
2869 +  mat - the matrix
2870 -  v - the location of the values
2871 
2872    Fortran Note:
2873    This routine is used differently from Fortran
2874 $    Mat         mat
2875 $    Scalar      mat_array(1)
2876 $    PetscOffset i_mat
2877 $    int         ierr
2878 $       call MatGetArray(mat,mat_array,i_mat,ierr)
2879 $
2880 $   Access first local entry in matrix with
2881 $      value = mat_array(i_mat + 1)
2882 $   (note here the array is treated as one dimensional)
2883 $
2884 $      ...... other code
2885 $       call MatRestoreArray(mat,mat_array,i_mat,ierr)
2886 
2887    See the Fortran chapter of the users manual and
2888    petsc/src/mat/examples/tests for details
2889 
2890 .keywords: matrix, array, elements, values, restore
2891 
2892 .seealso: MatGetArray()
2893 @*/
2894 int MatRestoreArray(Mat mat,Scalar **v)
2895 {
2896   int ierr;
2897 
2898   PetscFunctionBegin;
2899   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2900   PetscValidPointer(v);
2901   if (!mat->ops->restorearray) SETERRQ(PETSC_ERR_SUP,0,"");
2902   ierr = (*mat->ops->restorearray)(mat,v);CHKERRQ(ierr);
2903   PetscFunctionReturn(0);
2904 }
2905 
2906 #undef __FUNC__
2907 #define __FUNC__ "MatGetSubMatrices"
2908 /*@C
2909    MatGetSubMatrices - Extracts several submatrices from a matrix. If submat
2910    points to an array of valid matrices, they may be reused to store the new
2911    submatrices.
2912 
2913    Collective on Mat
2914 
2915    Input Parameters:
2916 +  mat - the matrix
2917 .  n   - the number of submatrixes to be extracted
2918 .  irow, icol - index sets of rows and columns to extract
2919 -  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
2920 
2921    Output Parameter:
2922 .  submat - the array of submatrices
2923 
2924    Notes:
2925    MatGetSubMatrices() can extract only sequential submatrices
2926    (from both sequential and parallel matrices). Use MatGetSubMatrix()
2927    to extract a parallel submatrix.
2928 
2929    When extracting submatrices from a parallel matrix, each processor can
2930    form a different submatrix by setting the rows and columns of its
2931    individual index sets according to the local submatrix desired.
2932 
2933    When finished using the submatrices, the user should destroy
2934    them with MatDestroySubMatrices().
2935 
2936    MAT_REUSE_MATRIX can only be used when the nonzero structure of the
2937    original matrix has not changed from that last call to MatGetSubMatrices()
2938 
2939    Fortran Note:
2940    The Fortran interface is slightly different from that given below, it
2941    requires one to pass in  as submat a Mat (integer) array of size at least m.
2942 
2943 .keywords: matrix, get, submatrix, submatrices
2944 
2945 .seealso: MatDestroyMatrices(), MatGetSubMatrix()
2946 @*/
2947 int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatReuse scall,Mat **submat)
2948 {
2949   int ierr;
2950 
2951   PetscFunctionBegin;
2952   PetscValidHeaderSpecific(mat,MAT_COOKIE);
2953   if (!mat->ops->getsubmatrices) SETERRQ(PETSC_ERR_SUP,0,"");
2954   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
2955 
2956   PLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);
2957   ierr = (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat); CHKERRQ(ierr);
2958   PLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);
2959 
2960   PetscFunctionReturn(0);
2961 }
2962 
2963 #undef __FUNC__
2964 #define __FUNC__ "MatDestroyMatrices"
2965 /*@C
2966    MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices().
2967 
2968    Collective on Mat
2969 
2970    Input Parameters:
2971 +  n - the number of local matrices
2972 -  mat - the matrices
2973 
2974 .keywords: matrix, destroy, submatrix, submatrices
2975 
2976 .seealso: MatGetSubMatrices()
2977 @*/
2978 int MatDestroyMatrices(int n,Mat **mat)
2979 {
2980   int ierr,i;
2981 
2982   PetscFunctionBegin;
2983   if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,1,"Trying to destroy negative number of matrices %d",n);
2984   PetscValidPointer(mat);
2985   for ( i=0; i<n; i++ ) {
2986     ierr = MatDestroy((*mat)[i]); CHKERRQ(ierr);
2987   }
2988   if (n) PetscFree(*mat);
2989   PetscFunctionReturn(0);
2990 }
2991 
2992 #undef __FUNC__
2993 #define __FUNC__ "MatIncreaseOverlap"
2994 /*@
2995    MatIncreaseOverlap - Given a set of submatrices indicated by index sets,
2996    replaces the index sets by larger ones that represent submatrices with
2997    additional overlap.
2998 
2999    Collective on Mat
3000 
3001    Input Parameters:
3002 +  mat - the matrix
3003 .  n   - the number of index sets
3004 .  is  - the array of pointers to index sets
3005 -  ov  - the additional overlap requested
3006 
3007 .keywords: matrix, overlap, Schwarz
3008 
3009 .seealso: MatGetSubMatrices()
3010 @*/
3011 int MatIncreaseOverlap(Mat mat,int n, IS *is,int ov)
3012 {
3013   int ierr;
3014 
3015   PetscFunctionBegin;
3016   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3017   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix");
3018   if (mat->factor)     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix");
3019 
3020   if (ov == 0) PetscFunctionReturn(0);
3021   if (!mat->ops->increaseoverlap) SETERRQ(PETSC_ERR_SUP,0,"");
3022   PLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);
3023   ierr = (*mat->ops->increaseoverlap)(mat,n,is,ov); CHKERRQ(ierr);
3024   PLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);
3025   PetscFunctionReturn(0);
3026 }
3027 
3028 #undef __FUNC__
3029 #define __FUNC__ "MatPrintHelp"
3030 /*@
3031    MatPrintHelp - Prints all the options for the matrix.
3032 
3033    Collective on Mat
3034 
3035    Input Parameter:
3036 .  mat - the matrix
3037 
3038    Options Database Keys:
3039 +  -help - Prints matrix options
3040 -  -h - Prints matrix options
3041 
3042 .keywords: mat, help
3043 
3044 .seealso: MatCreate(), MatCreateXXX()
3045 @*/
3046 int MatPrintHelp(Mat mat)
3047 {
3048   static int called = 0;
3049   int        ierr;
3050   MPI_Comm   comm;
3051 
3052   PetscFunctionBegin;
3053   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3054 
3055   comm = mat->comm;
3056   if (!called) {
3057     (*PetscHelpPrintf)(comm,"General matrix options:\n");
3058     (*PetscHelpPrintf)(comm,"  -mat_view_info: view basic matrix info during MatAssemblyEnd()\n");
3059     (*PetscHelpPrintf)(comm,"  -mat_view_info_detailed: view detailed matrix info during MatAssemblyEnd()\n");
3060     (*PetscHelpPrintf)(comm,"  -mat_view_draw: draw nonzero matrix structure during MatAssemblyEnd()\n");
3061     (*PetscHelpPrintf)(comm,"      -draw_pause <sec>: set seconds of display pause\n");
3062     (*PetscHelpPrintf)(comm,"      -display <name>: set alternate display\n");
3063     called = 1;
3064   }
3065   if (mat->ops->printhelp) {
3066     ierr = (*mat->ops->printhelp)(mat); CHKERRQ(ierr);
3067   }
3068   PetscFunctionReturn(0);
3069 }
3070 
3071 #undef __FUNC__
3072 #define __FUNC__ "MatGetBlockSize"
3073 /*@
3074    MatGetBlockSize - Returns the matrix block size; useful especially for the
3075    block row and block diagonal formats.
3076 
3077    Not Collective
3078 
3079    Input Parameter:
3080 .  mat - the matrix
3081 
3082    Output Parameter:
3083 .  bs - block size
3084 
3085    Notes:
3086    Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG.
3087    Block row formats are MATSEQBAIJ, MATMPIBAIJ
3088 
3089 .keywords: matrix, get, block, size
3090 
3091 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag()
3092 @*/
3093 int MatGetBlockSize(Mat mat,int *bs)
3094 {
3095   int ierr;
3096 
3097   PetscFunctionBegin;
3098   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3099   PetscValidIntPointer(bs);
3100   if (!mat->ops->getblocksize) SETERRQ(PETSC_ERR_SUP,0,"");
3101   ierr = (*mat->ops->getblocksize)(mat,bs);CHKERRQ(ierr);
3102   PetscFunctionReturn(0);
3103 }
3104 
3105 #undef __FUNC__
3106 #define __FUNC__ "MatGetRowIJ"
3107 /*@C
3108     MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices.
3109     EXPERTS ONLY.
3110 
3111    Collective on Mat
3112 
3113     Input Parameters:
3114 +   mat - the matrix
3115 .   shift -  0 or 1 indicating we want the indices starting at 0 or 1
3116 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3117                 symmetrized
3118 
3119     Output Parameters:
3120 +   n - number of rows in the (possibly compressed) matrix
3121 .   ia - the row pointers
3122 .   ja - the column indices
3123 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3124 
3125 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3126 @*/
3127 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3128 {
3129   int ierr;
3130 
3131   PetscFunctionBegin;
3132   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3133   if (ia) PetscValidIntPointer(ia);
3134   if (ja) PetscValidIntPointer(ja);
3135   PetscValidIntPointer(done);
3136   if (!mat->ops->getrowij) *done = PETSC_FALSE;
3137   else {
3138     *done = PETSC_TRUE;
3139     ierr  = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr);
3140   }
3141   PetscFunctionReturn(0);
3142 }
3143 
3144 #undef __FUNC__
3145 #define __FUNC__ "MatGetColumnIJ"
3146 /*@C
3147     MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices.
3148     EXPERTS ONLY.
3149 
3150     Collective on Mat
3151 
3152     Input Parameters:
3153 +   mat - the matrix
3154 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3155 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3156                 symmetrized
3157 
3158     Output Parameters:
3159 +   n - number of columns in the (possibly compressed) matrix
3160 .   ia - the column pointers
3161 .   ja - the row indices
3162 -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned
3163 
3164 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3165 @*/
3166 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3167 {
3168   int ierr;
3169 
3170   PetscFunctionBegin;
3171   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3172   if (ia) PetscValidIntPointer(ia);
3173   if (ja) PetscValidIntPointer(ja);
3174   PetscValidIntPointer(done);
3175 
3176   if (!mat->ops->getcolumnij) *done = PETSC_FALSE;
3177   else {
3178     *done = PETSC_TRUE;
3179     ierr  = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr);
3180   }
3181   PetscFunctionReturn(0);
3182 }
3183 
3184 #undef __FUNC__
3185 #define __FUNC__ "MatRestoreRowIJ"
3186 /*@C
3187     MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with
3188     MatGetRowIJ(). EXPERTS ONLY.
3189 
3190     Collective on Mat
3191 
3192     Input Parameters:
3193 +   mat - the matrix
3194 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3195 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3196                 symmetrized
3197 
3198     Output Parameters:
3199 +   n - size of (possibly compressed) matrix
3200 .   ia - the row pointers
3201 .   ja - the column indices
3202 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3203 
3204 .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
3205 @*/
3206 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3207 {
3208   int ierr;
3209 
3210   PetscFunctionBegin;
3211   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3212   if (ia) PetscValidIntPointer(ia);
3213   if (ja) PetscValidIntPointer(ja);
3214   PetscValidIntPointer(done);
3215 
3216   if (!mat->ops->restorerowij) *done = PETSC_FALSE;
3217   else {
3218     *done = PETSC_TRUE;
3219     ierr  = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr);
3220   }
3221   PetscFunctionReturn(0);
3222 }
3223 
3224 #undef __FUNC__
3225 #define __FUNC__ "MatRestoreColumnIJ"
3226 /*@C
3227     MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with
3228     MatGetColumnIJ(). EXPERTS ONLY.
3229 
3230     Collective on Mat
3231 
3232     Input Parameters:
3233 +   mat - the matrix
3234 .   shift - 1 or zero indicating we want the indices starting at 0 or 1
3235 -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
3236                 symmetrized
3237 
3238     Output Parameters:
3239 +   n - size of (possibly compressed) matrix
3240 .   ia - the column pointers
3241 .   ja - the row indices
3242 -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned
3243 
3244 .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
3245 @*/
3246 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
3247 {
3248   int ierr;
3249 
3250   PetscFunctionBegin;
3251   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3252   if (ia) PetscValidIntPointer(ia);
3253   if (ja) PetscValidIntPointer(ja);
3254   PetscValidIntPointer(done);
3255 
3256   if (!mat->ops->restorecolumnij) *done = PETSC_FALSE;
3257   else {
3258     *done = PETSC_TRUE;
3259     ierr  = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr);
3260   }
3261   PetscFunctionReturn(0);
3262 }
3263 
3264 #undef __FUNC__
3265 #define __FUNC__ "MatColoringPatch"
3266 /*@C
3267     MatColoringPatch - EXPERTS ONLY, used inside matrix coloring routines that
3268     use MatGetRowIJ() and/or MatGetColumnIJ().
3269 
3270     Collective on Mat
3271 
3272     Input Parameters:
3273 +   mat - the matrix
3274 .   n   - number of colors
3275 -   colorarray - array indicating color for each column
3276 
3277     Output Parameters:
3278 .   iscoloring - coloring generated using colorarray information
3279 
3280 .seealso: MatGetRowIJ(), MatGetColumnIJ()
3281 
3282 .keywords: mat, coloring, patch
3283 @*/
3284 int MatColoringPatch(Mat mat,int n,int *colorarray,ISColoring *iscoloring)
3285 {
3286   int ierr;
3287 
3288   PetscFunctionBegin;
3289   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3290   PetscValidIntPointer(colorarray);
3291 
3292   if (!mat->ops->coloringpatch) {SETERRQ(PETSC_ERR_SUP,0,"");}
3293   else {
3294     ierr  = (*mat->ops->coloringpatch)(mat,n,colorarray,iscoloring); CHKERRQ(ierr);
3295   }
3296   PetscFunctionReturn(0);
3297 }
3298 
3299 
3300 #undef __FUNC__
3301 #define __FUNC__ "MatSetUnfactored"
3302 /*@
3303    MatSetUnfactored - Resets a factored matrix to be treated as unfactored.
3304 
3305    Collective on Mat
3306 
3307    Input Parameter:
3308 .  mat - the factored matrix to be reset
3309 
3310    Notes:
3311    This routine should be used only with factored matrices formed by in-place
3312    factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE
3313    format).  This option can save memory, for example, when solving nonlinear
3314    systems with a matrix-free Newton-Krylov method and a matrix-based, in-place
3315    ILU(0) preconditioner.
3316 
3317   Note that one can specify in-place ILU(0) factorization by calling
3318 $     PCType(pc,PCILU);
3319 $     PCILUSeUseInPlace(pc);
3320   or by using the options -pc_type ilu -pc_ilu_in_place
3321 
3322   In-place factorization ILU(0) can also be used as a local
3323   solver for the blocks within the block Jacobi or additive Schwarz
3324   methods (runtime option: -sub_pc_ilu_in_place).  See the discussion
3325   of these preconditioners in the users manual for details on setting
3326   local solver options.
3327 
3328 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace()
3329 
3330 .keywords: matrix-free, in-place ILU, in-place LU
3331 @*/
3332 int MatSetUnfactored(Mat mat)
3333 {
3334   int ierr;
3335 
3336   PetscFunctionBegin;
3337   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3338   mat->factor = 0;
3339   if (!mat->ops->setunfactored) PetscFunctionReturn(0);
3340   ierr = (*mat->ops->setunfactored)(mat); CHKERRQ(ierr);
3341   PetscFunctionReturn(0);
3342 }
3343 
3344 #undef __FUNC__
3345 #define __FUNC__ "MatGetType"
3346 /*@C
3347    MatGetType - Gets the matrix type and name (as a string) from the matrix.
3348 
3349    Not Collective
3350 
3351    Input Parameter:
3352 .  mat - the matrix
3353 
3354    Output Parameter:
3355 +  type - the matrix type (or use PETSC_NULL)
3356 -  name - name of matrix type (or use PETSC_NULL)
3357 
3358 .keywords: matrix, get, type, name
3359 @*/
3360 int MatGetType(Mat mat,MatType *type,char **name)
3361 {
3362   int  itype = (int)mat->type;
3363   char *matname[10];
3364 
3365   PetscFunctionBegin;
3366   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3367 
3368   if (type) *type = (MatType) mat->type;
3369   if (name) {
3370     /* Note:  Be sure that this list corresponds to the enum in mat.h */
3371     matname[0] = "MATSEQDENSE";
3372     matname[1] = "MATSEQAIJ";
3373     matname[2] = "MATMPIAIJ";
3374     matname[3] = "MATSHELL";
3375     matname[4] = "MATMPIROWBS";
3376     matname[5] = "MATSEQBDIAG";
3377     matname[6] = "MATMPIBDIAG";
3378     matname[7] = "MATMPIDENSE";
3379     matname[8] = "MATSEQBAIJ";
3380     matname[9] = "MATMPIBAIJ";
3381 
3382     if (itype < 0 || itype > 9) *name = "Unknown matrix type";
3383     else                        *name = matname[itype];
3384   }
3385   PetscFunctionReturn(0);
3386 }
3387 
3388 /*MC
3389     MatGetArrayF90 - Accesses a matrix array from Fortran90.
3390 
3391     Synopsis:
3392     MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3393 
3394     Not collective
3395 
3396     Input Parameter:
3397 .   x - matrix
3398 
3399     Output Parameters:
3400 +   xx_v - the Fortran90 pointer to the array
3401 -   ierr - error code
3402 
3403     Example of Usage:
3404 .vb
3405       Scalar, pointer xx_v(:)
3406       ....
3407       call MatGetArrayF90(x,xx_v,ierr)
3408       a = xx_v(3)
3409       call MatRestoreArrayF90(x,xx_v,ierr)
3410 .ve
3411 
3412     Notes:
3413      Not yet supported for all F90 compilers
3414 
3415 .seealso:  MatRestoreArrayF90(), MatGetArray(), MatRestoreArray()
3416 
3417 .keywords:  matrix, array, f90
3418 M*/
3419 
3420 /*MC
3421     MatRestoreArrayF90 - Restores a matrix array that has been
3422     accessed with MatGetArrayF90().
3423 
3424     Synopsis:
3425     MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)
3426 
3427     Not collective
3428 
3429     Input Parameters:
3430 +   x - matrix
3431 -   xx_v - the Fortran90 pointer to the array
3432 
3433     Output Parameter:
3434 .   ierr - error code
3435 
3436     Example of Usage:
3437 .vb
3438        Scalar, pointer xx_v(:)
3439        ....
3440        call MatGetArrayF90(x,xx_v,ierr)
3441        a = xx_v(3)
3442        call MatRestoreArrayF90(x,xx_v,ierr)
3443 .ve
3444 
3445     Notes:
3446      Not yet supported for all F90 compilers
3447 
3448 .seealso:  MatGetArrayF90(), MatGetArray(), MatRestoreArray()
3449 
3450 .keywords:  matrix, array, f90
3451 M*/
3452 
3453 
3454 #undef __FUNC__
3455 #define __FUNC__ "MatGetSubMatrix"
3456 /*@
3457     MatGetSubMatrix - Gets a single submatrix on the same number of processors
3458                       as the original matrix.
3459 
3460    Collective on Mat
3461 
3462     Input Parameters:
3463 +   mat - the original matrix
3464 .   isrow - rows this processor should obtain
3465 .   iscol - columns for all processors you wish to keep
3466 .   csize - number of columns "local" to this processor (does nothing for sequential
3467             matrices). This should match the result from VecGetLocalSize(x,...) if you
3468             plan to use the matrix in a A*x or use PETSC_DECIDE
3469 -  cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
3470 
3471     Output Parameter:
3472 .   newmat - the new submatrix, of the same type as the old
3473 
3474 .seealso: MatGetSubMatrices()
3475 
3476 @*/
3477 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat)
3478 {
3479   int     ierr, size;
3480   Mat     *local;
3481 
3482   PetscFunctionBegin;
3483   MPI_Comm_size(mat->comm,&size);
3484 
3485   /* if original matrix is on just one processor then use submatrix generated */
3486   if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) {
3487     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr);
3488     PetscFunctionReturn(0);
3489   } else if (!mat->ops->getsubmatrix && size == 1) {
3490     ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr);
3491     *newmat = *local;
3492     PetscFree(local);
3493     PetscFunctionReturn(0);
3494   }
3495 
3496   if (!mat->ops->getsubmatrix) SETERRQ(PETSC_ERR_SUP,0,"Not currently implemented");
3497   ierr = (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr);
3498   PetscFunctionReturn(0);
3499 }
3500 
3501 #undef __FUNC__
3502 #define __FUNC__ "MatGetMaps"
3503 /*@C
3504    MatGetMaps - Returns the maps associated with the matrix.
3505 
3506    Not Collective
3507 
3508    Input Parameter:
3509 .  mat - the matrix
3510 
3511    Output Parameters:
3512 +  rmap - the row (right) map
3513 -  cmap - the column (left) map
3514 
3515 .keywords: matrix, get, map
3516 @*/
3517 int MatGetMaps(Mat mat,Map *rmap,Map *cmap)
3518 {
3519   int ierr;
3520 
3521   PetscFunctionBegin;
3522   PetscValidHeaderSpecific(mat,MAT_COOKIE);
3523   ierr = (*mat->ops->getmaps)(mat,rmap,cmap);CHKERRQ(ierr);
3524   PetscFunctionReturn(0);
3525 }
3526 
3527 /*
3528       Version that works for all PETSc matrices
3529 */
3530 #undef __FUNC__
3531 #define __FUNC__ "MatGetMaps_Petsc"
3532 int MatGetMaps_Petsc(Mat mat,Map *rmap,Map *cmap)
3533 {
3534   PetscFunctionBegin;
3535   if (rmap) *rmap = mat->rmap;
3536   if (cmap) *cmap = mat->cmap;
3537   PetscFunctionReturn(0);
3538 }
3539