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