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