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