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