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