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