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