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