Actual source code: matrix.c

  1: /*$Id: matrix.c,v 1.414 2001/09/28 18:57:28 balay Exp $*/

  3: /*
  4:    This is where the abstract matrix operations are defined
  5: */

 7:  #include src/mat/matimpl.h
 8:  #include src/vec/vecimpl.h

 10: /* Logging support */
 11: int MAT_COOKIE;
 12: int MATSNESMFCTX_COOKIE;
 13: int MAT_Mult, MAT_MultMatrixFree, MAT_Mults, MAT_MultConstrained, MAT_MultAdd, MAT_MultTranspose;
 14: int MAT_MultTransposeConstrained, MAT_MultTransposeAdd, MAT_Solve, MAT_Solves, MAT_SolveAdd, MAT_SolveTranspose;
 15: int MAT_SolveTransposeAdd, MAT_Relax, MAT_ForwardSolve, MAT_BackwardSolve, MAT_LUFactor, MAT_LUFactorSymbolic;
 16: int MAT_LUFactorNumeric, MAT_CholeskyFactor, MAT_CholeskyFactorSymbolic, MAT_CholeskyFactorNumeric, MAT_ILUFactor;
 17: int MAT_ILUFactorSymbolic, MAT_ICCFactorSymbolic, MAT_Copy, MAT_Convert, MAT_Scale, MAT_AssemblyBegin;
 18: int MAT_AssemblyEnd, MAT_SetValues, MAT_GetValues, MAT_GetRow, MAT_GetSubMatrices, MAT_GetColoring, MAT_GetOrdering;
 19: int MAT_IncreaseOverlap, MAT_Partitioning, MAT_ZeroEntries, MAT_Load, MAT_View, MAT_AXPY, MAT_FDColoringCreate;
 20: int MAT_FDColoringApply,MAT_Transpose;

 22: #undef __FUNCT__  
 24: /*@C
 25:    MatGetRow - Gets a row of a matrix.  You MUST call MatRestoreRow()
 26:    for each row that you get to ensure that your application does
 27:    not bleed memory.

 29:    Not Collective

 31:    Input Parameters:
 32: +  mat - the matrix
 33: -  row - the row to get

 35:    Output Parameters:
 36: +  ncols -  the number of nonzeros in the row
 37: .  cols - if not NULL, the column numbers
 38: -  vals - if not NULL, the values

 40:    Notes:
 41:    This routine is provided for people who need to have direct access
 42:    to the structure of a matrix.  We hope that we provide enough
 43:    high-level matrix routines that few users will need it. 

 45:    MatGetRow() always returns 0-based column indices, regardless of
 46:    whether the internal representation is 0-based (default) or 1-based.

 48:    For better efficiency, set cols and/or vals to PETSC_NULL if you do
 49:    not wish to extract these quantities.

 51:    The user can only examine the values extracted with MatGetRow();
 52:    the values cannot be altered.  To change the matrix entries, one
 53:    must use MatSetValues().

 55:    You can only have one call to MatGetRow() outstanding for a particular
 56:    matrix at a time, per processor. MatGetRow() can only obtained rows
 57:    associated with the given processor, it cannot get rows from the 
 58:    other processors; for that we suggest using MatGetSubMatrices(), then
 59:    MatGetRow() on the submatrix. The row indix passed to MatGetRows() 
 60:    is in the global number of rows.

 62:    Fortran Notes:
 63:    The calling sequence from Fortran is 
 64: .vb
 65:    MatGetRow(matrix,row,ncols,cols,values,ierr)
 66:          Mat     matrix (input)
 67:          integer row    (input)
 68:          integer ncols  (output)
 69:          integer cols(maxcols) (output)
 70:          double precision (or double complex) values(maxcols) output
 71: .ve
 72:    where maxcols >= maximum nonzeros in any row of the matrix.

 74:    Caution:
 75:    Do not try to change the contents of the output arrays (cols and vals).
 76:    In some cases, this may corrupt the matrix.

 78:    Level: advanced

 80:    Concepts: matrices^row access

 82: .seealso: MatRestoreRow(), MatSetValues(), MatGetValues(), MatGetSubmatrices(), MatGetDiagonal()
 83: @*/
 84: int MatGetRow(Mat mat,int row,int *ncols,int **cols,PetscScalar **vals)
 85: {
 86:   int   ierr;

 92:   MatPreallocated(mat);
 93:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
 94:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
 95:   if (!mat->ops->getrow) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
 96:   PetscLogEventBegin(MAT_GetRow,mat,0,0,0);
 97:   (*mat->ops->getrow)(mat,row,ncols,cols,vals);
 98:   PetscLogEventEnd(MAT_GetRow,mat,0,0,0);
 99:   return(0);
100: }

102: #undef __FUNCT__  
104: /*@C  
105:    MatRestoreRow - Frees any temporary space allocated by MatGetRow().

107:    Not Collective

109:    Input Parameters:
110: +  mat - the matrix
111: .  row - the row to get
112: .  ncols, cols - the number of nonzeros and their columns
113: -  vals - if nonzero the column values

115:    Notes: 
116:    This routine should be called after you have finished examining the entries.

118:    Fortran Notes:
119:    The calling sequence from Fortran is 
120: .vb
121:    MatRestoreRow(matrix,row,ncols,cols,values,ierr)
122:       Mat     matrix (input)
123:       integer row    (input)
124:       integer ncols  (output)
125:       integer cols(maxcols) (output)
126:       double precision (or double complex) values(maxcols) output
127: .ve
128:    Where maxcols >= maximum nonzeros in any row of the matrix. 

130:    In Fortran MatRestoreRow() MUST be called after MatGetRow()
131:    before another call to MatGetRow() can be made.

133:    Level: advanced

135: .seealso:  MatGetRow()
136: @*/
137: int MatRestoreRow(Mat mat,int row,int *ncols,int **cols,PetscScalar **vals)
138: {

144:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
145:   if (!mat->ops->restorerow) return(0);
146:   (*mat->ops->restorerow)(mat,row,ncols,cols,vals);
147:   return(0);
148: }

150: #undef __FUNCT__  
152: /*@C
153:    MatView - Visualizes a matrix object.

155:    Collective on Mat

157:    Input Parameters:
158: +  mat - the matrix
159: -  ptr - visualization context

161:   Notes:
162:   The available visualization contexts include
163: +    PETSC_VIEWER_STDOUT_SELF - standard output (default)
164: .    PETSC_VIEWER_STDOUT_WORLD - synchronized standard
165:         output where only the first processor opens
166:         the file.  All other processors send their 
167:         data to the first processor to print. 
168: -     PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure

170:    The user can open alternative visualization contexts with
171: +    PetscViewerASCIIOpen() - Outputs matrix to a specified file
172: .    PetscViewerBinaryOpen() - Outputs matrix in binary to a
173:          specified file; corresponding input uses MatLoad()
174: .    PetscViewerDrawOpen() - Outputs nonzero matrix structure to 
175:          an X window display
176: -    PetscViewerSocketOpen() - Outputs matrix to Socket viewer.
177:          Currently only the sequential dense and AIJ
178:          matrix types support the Socket viewer.

180:    The user can call PetscViewerSetFormat() to specify the output
181:    format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
182:    PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen).  Available formats include
183: +    PETSC_VIEWER_ASCII_DEFAULT - default, prints matrix contents
184: .    PETSC_VIEWER_ASCII_MATLAB - prints matrix contents in Matlab format
185: .    PETSC_VIEWER_ASCII_DENSE - prints entire matrix including zeros
186: .    PETSC_VIEWER_ASCII_COMMON - prints matrix contents, using a sparse 
187:          format common among all matrix types
188: .    PETSC_VIEWER_ASCII_IMPL - prints matrix contents, using an implementation-specific 
189:          format (which is in many cases the same as the default)
190: .    PETSC_VIEWER_ASCII_INFO - prints basic information about the matrix
191:          size and structure (not the matrix entries)
192: .    PETSC_VIEWER_ASCII_INFO_DETAIL - prints more detailed information about
193:          the matrix structure

195:    Options Database Keys:
196: +  -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
197: .  -mat_view_info_detailed - Prints more detailed info
198: .  -mat_view - Prints matrix in ASCII format
199: .  -mat_view_matlab - Prints matrix in Matlab format
200: .  -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
201: .  -display <name> - Sets display name (default is host)
202: .  -draw_pause <sec> - Sets number of seconds to pause after display
203: .  -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
204: .  -viewer_socket_machine <machine>
205: .  -viewer_socket_port <port>
206: .  -mat_view_binary - save matrix to file in binary format
207: -  -viewer_binary_filename <name>
208:    Level: beginner

210:    Concepts: matrices^viewing
211:    Concepts: matrices^plotting
212:    Concepts: matrices^printing

214: .seealso: PetscViewerSetFormat(), PetscViewerASCIIOpen(), PetscViewerDrawOpen(), 
215:           PetscViewerSocketOpen(), PetscViewerBinaryOpen(), MatLoad()
216: @*/
217: int MatView(Mat mat,PetscViewer viewer)
218: {
219:   int               ierr,rows,cols;
220:   PetscTruth        isascii;
221:   char              *cstr;
222:   PetscViewerFormat format;

227:   MatPreallocated(mat);
228:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(mat->comm);
231:   if (!mat->assembled) SETERRQ(1,"Must call MatAssemblyBegin/End() before viewing matrix");

233:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
234:   if (isascii) {
235:     PetscViewerGetFormat(viewer,&format);
236:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
237:       if (mat->prefix) {
238:         PetscViewerASCIIPrintf(viewer,"Matrix Object:(%s)n",mat->prefix);
239:       } else {
240:         PetscViewerASCIIPrintf(viewer,"Matrix Object:n");
241:       }
242:       PetscViewerASCIIPushTab(viewer);
243:       MatGetType(mat,&cstr);
244:       MatGetSize(mat,&rows,&cols);
245:       PetscViewerASCIIPrintf(viewer,"type=%s, rows=%d, cols=%dn",cstr,rows,cols);
246:       if (mat->ops->getinfo) {
247:         MatInfo info;
248:         MatGetInfo(mat,MAT_GLOBAL_SUM,&info);
249:         PetscViewerASCIIPrintf(viewer,"total: nonzeros=%d, allocated nonzeros=%dn",
250:                           (int)info.nz_used,(int)info.nz_allocated);
251:       }
252:     }
253:   }
254:   if (mat->ops->view) {
255:     PetscViewerASCIIPushTab(viewer);
256:     (*mat->ops->view)(mat,viewer);
257:     PetscViewerASCIIPopTab(viewer);
258:   } else if (!isascii) {
259:     SETERRQ1(1,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
260:   }
261:   if (isascii) {
262:     PetscViewerGetFormat(viewer,&format);
263:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
264:       PetscViewerASCIIPopTab(viewer);
265:     }
266:   }
267:   return(0);
268: }

270: #undef __FUNCT__  
272: /*@C
273:    MatScaleSystem - Scale a vector solution and right hand side to 
274:    match the scaling of a scaled matrix.
275:   
276:    Collective on Mat

278:    Input Parameter:
279: +  mat - the matrix
280: .  x - solution vector (or PETSC_NULL)
281: -  b - right hand side vector (or PETSC_NULL)


284:    Notes: 
285:    For AIJ, BAIJ, and BDiag matrix formats, the matrices are not 
286:    internally scaled, so this does nothing. For MPIROWBS it
287:    permutes and diagonally scales.

289:    The SLES methods automatically call this routine when required
290:    (via PCPreSolve()) so it is rarely used directly.

292:    Level: Developer            

294:    Concepts: matrices^scaling

296: .seealso: MatUseScaledForm(), MatUnScaleSystem()
297: @*/
298: int MatScaleSystem(Mat mat,Vec x,Vec b)
299: {

305:   MatPreallocated(mat);

309:   if (mat->ops->scalesystem) {
310:     (*mat->ops->scalesystem)(mat,x,b);
311:   }
312:   return(0);
313: }

315: #undef __FUNCT__  
317: /*@C
318:    MatUnScaleSystem - Unscales a vector solution and right hand side to 
319:    match the original scaling of a scaled matrix.
320:   
321:    Collective on Mat

323:    Input Parameter:
324: +  mat - the matrix
325: .  x - solution vector (or PETSC_NULL)
326: -  b - right hand side vector (or PETSC_NULL)


329:    Notes: 
330:    For AIJ, BAIJ, and BDiag matrix formats, the matrices are not 
331:    internally scaled, so this does nothing. For MPIROWBS it
332:    permutes and diagonally scales.

334:    The SLES methods automatically call this routine when required
335:    (via PCPreSolve()) so it is rarely used directly.

337:    Level: Developer            

339: .seealso: MatUseScaledForm(), MatScaleSystem()
340: @*/
341: int MatUnScaleSystem(Mat mat,Vec x,Vec b)
342: {

348:   MatPreallocated(mat);
351:   if (mat->ops->unscalesystem) {
352:     (*mat->ops->unscalesystem)(mat,x,b);
353:   }
354:   return(0);
355: }

357: #undef __FUNCT__  
359: /*@C
360:    MatUseScaledForm - For matrix storage formats that scale the 
361:    matrix (for example MPIRowBS matrices are diagonally scaled on
362:    assembly) indicates matrix operations (MatMult() etc) are 
363:    applied using the scaled matrix.
364:   
365:    Collective on Mat

367:    Input Parameter:
368: +  mat - the matrix
369: -  scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for 
370:             applying the original matrix

372:    Notes: 
373:    For scaled matrix formats, applying the original, unscaled matrix
374:    will be slightly more expensive

376:    Level: Developer            

378: .seealso: MatScaleSystem(), MatUnScaleSystem()
379: @*/
380: int MatUseScaledForm(Mat mat,PetscTruth scaled)
381: {

387:   MatPreallocated(mat);
388:   if (mat->ops->usescaledform) {
389:     (*mat->ops->usescaledform)(mat,scaled);
390:   }
391:   return(0);
392: }

394: #undef __FUNCT__  
396: /*@C
397:    MatDestroy - Frees space taken by a matrix.
398:   
399:    Collective on Mat

401:    Input Parameter:
402: .  A - the matrix

404:    Level: beginner

406: @*/
407: int MatDestroy(Mat A)
408: {

414:   MatPreallocated(A);
415:   if (--A->refct > 0) return(0);

417:   /* if memory was published with AMS then destroy it */
418:   PetscObjectDepublish(A);
419:   if (A->mapping) {
420:     ISLocalToGlobalMappingDestroy(A->mapping);
421:   }
422:   if (A->bmapping) {
423:     ISLocalToGlobalMappingDestroy(A->bmapping);
424:   }
425:   if (A->rmap) {
426:     PetscMapDestroy(A->rmap);
427:   }
428:   if (A->cmap) {
429:     PetscMapDestroy(A->cmap);
430:   }

432:   (*A->ops->destroy)(A);
433:   PetscLogObjectDestroy(A);
434:   PetscHeaderDestroy(A);
435:   return(0);
436: }

438: #undef __FUNCT__  
440: /*@
441:    MatValid - Checks whether a matrix object is valid.

443:    Collective on Mat

445:    Input Parameter:
446: .  m - the matrix to check 

448:    Output Parameter:
449:    flg - flag indicating matrix status, either
450:    PETSC_TRUE if matrix is valid, or PETSC_FALSE otherwise.

452:    Level: developer

454:    Concepts: matrices^validity
455: @*/
456: int MatValid(Mat m,PetscTruth *flg)
457: {
460:   if (!m)                           *flg = PETSC_FALSE;
461:   else if (m->cookie != MAT_COOKIE) *flg = PETSC_FALSE;
462:   else                              *flg = PETSC_TRUE;
463:   return(0);
464: }

466: #undef __FUNCT__  
468: /*@ 
469:    MatSetValues - Inserts or adds a block of values into a matrix.
470:    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 
471:    MUST be called after all calls to MatSetValues() have been completed.

473:    Not Collective

475:    Input Parameters:
476: +  mat - the matrix
477: .  v - a logically two-dimensional array of values
478: .  m, idxm - the number of rows and their global indices 
479: .  n, idxn - the number of columns and their global indices
480: -  addv - either ADD_VALUES or INSERT_VALUES, where
481:    ADD_VALUES adds values to any existing entries, and
482:    INSERT_VALUES replaces existing entries with new values

484:    Notes:
485:    By default the values, v, are row-oriented and unsorted.
486:    See MatSetOption() for other options.

488:    Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES 
489:    options cannot be mixed without intervening calls to the assembly
490:    routines.

492:    MatSetValues() uses 0-based row and column numbers in Fortran 
493:    as well as in C.

495:    Negative indices may be passed in idxm and idxn, these rows and columns are 
496:    simply ignored. This allows easily inserting element stiffness matrices
497:    with homogeneous Dirchlet boundary conditions that you don't want represented
498:    in the matrix.

500:    Efficiency Alert:
501:    The routine MatSetValuesBlocked() may offer much better efficiency
502:    for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).

504:    Level: beginner

506:    Concepts: matrices^putting entries in

508: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
509: @*/
510: int MatSetValues(Mat mat,int m,int *idxm,int n,int *idxn,PetscScalar *v,InsertMode addv)
511: {

515:   if (!m || !n) return(0); /* no values to insert */
518:   MatPreallocated(mat);
522:   if (mat->insertmode == NOT_SET_VALUES) {
523:     mat->insertmode = addv;
524:   }
525: #if defined(PETSC_USE_BOPT_g)
526:   else if (mat->insertmode != addv) {
527:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
528:   }
529:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
530: #endif

532:   if (mat->assembled) {
533:     mat->was_assembled = PETSC_TRUE;
534:     mat->assembled     = PETSC_FALSE;
535:   }
536:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
537:   if (!mat->ops->setvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
538:   (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);
539:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
540:   return(0);
541: }

543: #undef __FUNCT__  
545: /*@C 
546:    MatSetValuesStencil - Inserts or adds a block of values into a matrix.
547:      Using structured grid indexing

549:    Not Collective

551:    Input Parameters:
552: +  mat - the matrix
553: .  v - a logically two-dimensional array of values
554: .  m - number of rows being entered
555: .  idxm - grid coordinates (and component number when dof > 1) for matrix rows being entered
556: .  n - number of columns being entered
557: .  idxn - grid coordinates (and component number when dof > 1) for matrix columns being entered 
558: -  addv - either ADD_VALUES or INSERT_VALUES, where
559:    ADD_VALUES adds values to any existing entries, and
560:    INSERT_VALUES replaces existing entries with new values

562:    Notes:
563:    By default the values, v, are row-oriented and unsorted.
564:    See MatSetOption() for other options.

566:    Calls to MatSetValuesStencil() with the INSERT_VALUES and ADD_VALUES 
567:    options cannot be mixed without intervening calls to the assembly
568:    routines.

570:    The grid coordinates are across the entire grid, not just the local portion

572:    MatSetValuesStencil() uses 0-based row and column numbers in Fortran 
573:    as well as in C.

575:    In order to use this routine you must either obtain the matrix with DAGetMatrix()
576:    or call MatSetLocalToGlobalMapping() and MatSetStencil() first.

578:    In Fortran idxm and idxn should be declared as
579: $     MatStencil idxm(4,m),idxn(4,n)
580:    and the values inserted using
581: $    idxm(MatStencil_i,1) = i
582: $    idxm(MatStencil_j,1) = j
583: $    idxm(MatStencil_k,1) = k
584: $    idxm(MatStencil_c,1) = c
585:    etc

587:    Negative indices may be passed in idxm and idxn, these rows and columns are 
588:    simply ignored. This allows easily inserting element stiffness matrices
589:    with homogeneous Dirchlet boundary conditions that you don't want represented
590:    in the matrix.

592:    Inspired by the structured grid interface to the HYPRE package
593:    (http://www.llnl.gov/CASC/hypre)

595:    Efficiency Alert:
596:    The routine MatSetValuesBlockedStencil() may offer much better efficiency
597:    for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).

599:    Level: beginner

601:    Concepts: matrices^putting entries in

603: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
604:           MatSetValues(), MatSetValuesBlockedStencil(), MatSetStencil(), DAGetMatrix()
605: @*/
606: int MatSetValuesStencil(Mat mat,int m,MatStencil *idxm,int n,MatStencil *idxn,PetscScalar *v,InsertMode addv)
607: {
608:   int j,i,ierr,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
609:   int *starts = mat->stencil.starts,*dxm = (int*)idxm,*dxn = (int*)idxn,sdim = dim - (1 - (int)mat->stencil.noc);

612:   if (!m || !n) return(0); /* no values to insert */

619:   if (m > 128) SETERRQ1(1,"Can only set 128 rows at a time; trying to set %d",m);
620:   if (n > 128) SETERRQ1(1,"Can only set 256 columns at a time; trying to set %d",n);

622:   for (i=0; i<m; i++) {
623:     for (j=0; j<3-sdim; j++) dxm++;
624:     tmp = *dxm++ - starts[0];
625:     for (j=0; j<dim-1; j++) {
626:       tmp = tmp*dims[j] + *dxm++ - starts[j+1];
627:     }
628:     if (mat->stencil.noc) dxm++;
629:     jdxm[i] = tmp;
630:   }
631:   for (i=0; i<n; i++) {
632:     for (j=0; j<3-sdim; j++) dxn++;
633:     tmp = *dxn++ - starts[0];
634:     for (j=0; j<dim-1; j++) {
635:       tmp = tmp*dims[j] + *dxn++ - starts[j+1];
636:     }
637:     if (mat->stencil.noc) dxn++;
638:     jdxn[i] = tmp;
639:   }
640:   MatSetValuesLocal(mat,m,jdxm,n,jdxn,v,addv);
641:   return(0);
642: }

644: #undef __FUNCT__  
646: /*@ 
647:    MatSetStencil - Sets the grid information for setting values into a matrix via
648:         MatSetStencil()

650:    Not Collective

652:    Input Parameters:
653: +  mat - the matrix
654: .  dim - dimension of the grid 1,2, or 3
655: .  dims - number of grid points in x, y, and z direction, including ghost points on your processor
656: .  starts - starting point of ghost nodes on your processor in x, y, and z direction 
657: -  dof - number of degrees of freedom per node


660:    Inspired by the structured grid interface to the HYPRE package
661:    (www.llnl.gov/CASC/hyper)

663:    Level: beginner

665:    Concepts: matrices^putting entries in

667: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
668:           MatSetValues(), MatSetValuesBlockedStencil(), MatSetValuesStencil()
669: @*/
670: int MatSetStencil(Mat mat,int dim,int *dims,int *starts,int dof)
671: {
672:   int i;


679:   mat->stencil.dim = dim + (dof > 1);
680:   for (i=0; i<dim; i++) {
681:     mat->stencil.dims[i]   = dims[dim-i-1];      /* copy the values in backwards */
682:     mat->stencil.starts[i] = starts[dim-i-1];
683:   }
684:   mat->stencil.dims[dim]   = dof;
685:   mat->stencil.starts[dim] = 0;
686:   mat->stencil.noc         = (PetscTruth)(dof == 1);
687:   return(0);
688: }

690: #undef __FUNCT__  
692: /*@ 
693:    MatSetValuesBlocked - Inserts or adds a block of values into a matrix.

695:    Not Collective

697:    Input Parameters:
698: +  mat - the matrix
699: .  v - a logically two-dimensional array of values
700: .  m, idxm - the number of block rows and their global block indices 
701: .  n, idxn - the number of block columns and their global block indices
702: -  addv - either ADD_VALUES or INSERT_VALUES, where
703:    ADD_VALUES adds values to any existing entries, and
704:    INSERT_VALUES replaces existing entries with new values

706:    Notes:
707:    By default the values, v, are row-oriented and unsorted. So the layout of 
708:    v is the same as for MatSetValues(). See MatSetOption() for other options.

710:    Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES 
711:    options cannot be mixed without intervening calls to the assembly
712:    routines.

714:    MatSetValuesBlocked() uses 0-based row and column numbers in Fortran 
715:    as well as in C.

717:    Negative indices may be passed in idxm and idxn, these rows and columns are 
718:    simply ignored. This allows easily inserting element stiffness matrices
719:    with homogeneous Dirchlet boundary conditions that you don't want represented
720:    in the matrix.

722:    Each time an entry is set within a sparse matrix via MatSetValues(),
723:    internal searching must be done to determine where to place the the
724:    data in the matrix storage space.  By instead inserting blocks of 
725:    entries via MatSetValuesBlocked(), the overhead of matrix assembly is
726:    reduced.

728:    Restrictions:
729:    MatSetValuesBlocked() is currently supported only for the block AIJ
730:    matrix format (MATSEQBAIJ and MATMPIBAIJ, which are created via
731:    MatCreateSeqBAIJ() and MatCreateMPIBAIJ()).

733:    Level: intermediate

735:    Concepts: matrices^putting entries in blocked

737: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
738: @*/
739: int MatSetValuesBlocked(Mat mat,int m,int *idxm,int n,int *idxn,PetscScalar *v,InsertMode addv)
740: {

744:   if (!m || !n) return(0); /* no values to insert */
747:   MatPreallocated(mat);
751:   if (mat->insertmode == NOT_SET_VALUES) {
752:     mat->insertmode = addv;
753:   }
754: #if defined(PETSC_USE_BOPT_g) 
755:   else if (mat->insertmode != addv) {
756:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
757:   }
758:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
759: #endif

761:   if (mat->assembled) {
762:     mat->was_assembled = PETSC_TRUE;
763:     mat->assembled     = PETSC_FALSE;
764:   }
765:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
766:   if (!mat->ops->setvaluesblocked) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
767:   (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);
768:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
769:   return(0);
770: }

772: /*MC
773:    MatSetValue - Set a single entry into a matrix.

775:    Synopsis:
776:    int MatSetValue(Mat m,int row,int col,PetscScalar value,InsertMode mode);

778:    Not collective

780:    Input Parameters:
781: +  m - the matrix
782: .  row - the row location of the entry
783: .  col - the column location of the entry
784: .  value - the value to insert
785: -  mode - either INSERT_VALUES or ADD_VALUES

787:    Notes: 
788:    For efficiency one should use MatSetValues() and set several or many
789:    values simultaneously if possible.

791:    Note that VecSetValue() does NOT return an error code (since this
792:    is checked internally).

794:    Level: beginner

796: .seealso: MatSetValues()
797: M*/

799: #undef __FUNCT__  
801: /*@ 
802:    MatGetValues - Gets a block of values from a matrix.

804:    Not Collective; currently only returns a local block

806:    Input Parameters:
807: +  mat - the matrix
808: .  v - a logically two-dimensional array for storing the values
809: .  m, idxm - the number of rows and their global indices 
810: -  n, idxn - the number of columns and their global indices

812:    Notes:
813:    The user must allocate space (m*n PetscScalars) for the values, v.
814:    The values, v, are then returned in a row-oriented format, 
815:    analogous to that used by default in MatSetValues().

817:    MatGetValues() uses 0-based row and column numbers in
818:    Fortran as well as in C.

820:    MatGetValues() requires that the matrix has been assembled
821:    with MatAssemblyBegin()/MatAssemblyEnd().  Thus, calls to
822:    MatSetValues() and MatGetValues() CANNOT be made in succession
823:    without intermediate matrix assembly.

825:    Level: advanced

827:    Concepts: matrices^accessing values

829: .seealso: MatGetRow(), MatGetSubmatrices(), MatSetValues()
830: @*/
831: int MatGetValues(Mat mat,int m,int *idxm,int n,int *idxn,PetscScalar *v)
832: {

838:   MatPreallocated(mat);
842:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
843:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
844:   if (!mat->ops->getvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

846:   PetscLogEventBegin(MAT_GetValues,mat,0,0,0);
847:   (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);
848:   PetscLogEventEnd(MAT_GetValues,mat,0,0,0);
849:   return(0);
850: }

852: #undef __FUNCT__  
854: /*@
855:    MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
856:    the routine MatSetValuesLocal() to allow users to insert matrix entries
857:    using a local (per-processor) numbering.

859:    Not Collective

861:    Input Parameters:
862: +  x - the matrix
863: -  mapping - mapping created with ISLocalToGlobalMappingCreate() 
864:              or ISLocalToGlobalMappingCreateIS()

866:    Level: intermediate

868:    Concepts: matrices^local to global mapping
869:    Concepts: local to global mapping^for matrices

871: .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
872: @*/
873: int MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping)
874: {
879:   MatPreallocated(x);
881:   if (x->mapping) {
882:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
883:   }

885:   if (x->ops->setlocaltoglobalmapping) {
886:     (*x->ops->setlocaltoglobalmapping)(x,mapping);
887:   } else {
888:     x->mapping = mapping;
889:     PetscObjectReference((PetscObject)mapping);
890:   }
891:   return(0);
892: }

894: #undef __FUNCT__  
896: /*@
897:    MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use
898:    by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
899:    entries using a local (per-processor) numbering.

901:    Not Collective

903:    Input Parameters:
904: +  x - the matrix
905: -  mapping - mapping created with ISLocalToGlobalMappingCreate() or
906:              ISLocalToGlobalMappingCreateIS()

908:    Level: intermediate

910:    Concepts: matrices^local to global mapping blocked
911:    Concepts: local to global mapping^for matrices, blocked

913: .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
914:            MatSetValuesBlocked(), MatSetValuesLocal()
915: @*/
916: int MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping mapping)
917: {
922:   MatPreallocated(x);
924:   if (x->bmapping) {
925:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
926:   }
927: 
928:   x->bmapping = mapping;
929:   PetscObjectReference((PetscObject)mapping);
930:   return(0);
931: }

933: #undef __FUNCT__  
935: /*@
936:    MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
937:    using a local ordering of the nodes. 

939:    Not Collective

941:    Input Parameters:
942: +  x - the matrix
943: .  nrow, irow - number of rows and their local indices
944: .  ncol, icol - number of columns and their local indices
945: .  y -  a logically two-dimensional array of values
946: -  addv - either INSERT_VALUES or ADD_VALUES, where
947:    ADD_VALUES adds values to any existing entries, and
948:    INSERT_VALUES replaces existing entries with new values

950:    Notes:
951:    Before calling MatSetValuesLocal(), the user must first set the
952:    local-to-global mapping by calling MatSetLocalToGlobalMapping().

954:    Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES 
955:    options cannot be mixed without intervening calls to the assembly
956:    routines.

958:    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 
959:    MUST be called after all calls to MatSetValuesLocal() have been completed.

961:    Level: intermediate

963:    Concepts: matrices^putting entries in with local numbering

965: .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping(),
966:            MatSetValueLocal()
967: @*/
968: int MatSetValuesLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,PetscScalar *y,InsertMode addv)
969: {
970:   int ierr,irowm[2048],icolm[2048];

975:   MatPreallocated(mat);

980:   if (mat->insertmode == NOT_SET_VALUES) {
981:     mat->insertmode = addv;
982:   }
983: #if defined(PETSC_USE_BOPT_g) 
984:   else if (mat->insertmode != addv) {
985:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
986:   }
987:   if (!mat->ops->setvalueslocal && (nrow > 2048 || ncol > 2048)) {
988:     SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
989:   }
990:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
991: #endif

993:   if (mat->assembled) {
994:     mat->was_assembled = PETSC_TRUE;
995:     mat->assembled     = PETSC_FALSE;
996:   }
997:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
998:   if (!mat->ops->setvalueslocal) {
999:     ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm);
1000:     ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);
1001:     (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);
1002:   } else {
1003:     (*mat->ops->setvalueslocal)(mat,nrow,irow,ncol,icol,y,addv);
1004:   }
1005:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1006:   return(0);
1007: }

1009: #undef __FUNCT__  
1011: /*@
1012:    MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
1013:    using a local ordering of the nodes a block at a time. 

1015:    Not Collective

1017:    Input Parameters:
1018: +  x - the matrix
1019: .  nrow, irow - number of rows and their local indices
1020: .  ncol, icol - number of columns and their local indices
1021: .  y -  a logically two-dimensional array of values
1022: -  addv - either INSERT_VALUES or ADD_VALUES, where
1023:    ADD_VALUES adds values to any existing entries, and
1024:    INSERT_VALUES replaces existing entries with new values

1026:    Notes:
1027:    Before calling MatSetValuesBlockedLocal(), the user must first set the
1028:    local-to-global mapping by calling MatSetLocalToGlobalMappingBlock(),
1029:    where the mapping MUST be set for matrix blocks, not for matrix elements.

1031:    Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES 
1032:    options cannot be mixed without intervening calls to the assembly
1033:    routines.

1035:    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 
1036:    MUST be called after all calls to MatSetValuesBlockedLocal() have been completed.

1038:    Level: intermediate

1040:    Concepts: matrices^putting blocked values in with local numbering

1042: .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked()
1043: @*/
1044: int MatSetValuesBlockedLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,PetscScalar *y,InsertMode addv)
1045: {
1046:   int ierr,irowm[2048],icolm[2048];

1051:   MatPreallocated(mat);
1055:   if (mat->insertmode == NOT_SET_VALUES) {
1056:     mat->insertmode = addv;
1057:   }
1058: #if defined(PETSC_USE_BOPT_g) 
1059:   else if (mat->insertmode != addv) {
1060:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
1061:   }
1062:   if (!mat->bmapping) {
1063:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with MatSetLocalToGlobalMappingBlock()");
1064:   }
1065:   if (nrow > 2048 || ncol > 2048) {
1066:     SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol);
1067:   }
1068:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1069: #endif

1071:   if (mat->assembled) {
1072:     mat->was_assembled = PETSC_TRUE;
1073:     mat->assembled     = PETSC_FALSE;
1074:   }
1075:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
1076:   ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm);
1077:   ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm);
1078:   (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);
1079:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
1080:   return(0);
1081: }

1083: /* --------------------------------------------------------*/
1084: #undef __FUNCT__  
1086: /*@
1087:    MatMult - Computes the matrix-vector product, y = Ax.

1089:    Collective on Mat and Vec

1091:    Input Parameters:
1092: +  mat - the matrix
1093: -  x   - the vector to be multilplied

1095:    Output Parameters:
1096: .  y - the result

1098:    Notes:
1099:    The vectors x and y cannot be the same.  I.e., one cannot
1100:    call MatMult(A,y,y).

1102:    Level: beginner

1104:    Concepts: matrix-vector product

1106: .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
1107: @*/
1108: int MatMult(Mat mat,Vec x,Vec y)
1109: {

1115:   MatPreallocated(mat);

1119:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1120:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1121:   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1122: #ifndef PETSC_HAVE_CONSTRAINTS
1123:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1124:   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
1125:   if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %d %d",mat->m,y->n);
1126: #endif

1128:   if (mat->nullsp) {
1129:     MatNullSpaceRemove(mat->nullsp,x,&x);
1130:   }

1132:   PetscLogEventBegin(MAT_Mult,mat,x,y,0);
1133:   (*mat->ops->mult)(mat,x,y);
1134:   PetscLogEventEnd(MAT_Mult,mat,x,y,0);

1136:   if (mat->nullsp) {
1137:     MatNullSpaceRemove(mat->nullsp,y,PETSC_NULL);
1138:   }
1139:   return(0);
1140: }

1142: #undef __FUNCT__  
1144: /*@
1145:    MatMultTranspose - Computes matrix transpose times a vector.

1147:    Collective on Mat and Vec

1149:    Input Parameters:
1150: +  mat - the matrix
1151: -  x   - the vector to be multilplied

1153:    Output Parameters:
1154: .  y - the result

1156:    Notes:
1157:    The vectors x and y cannot be the same.  I.e., one cannot
1158:    call MatMultTranspose(A,y,y).

1160:    Level: beginner

1162:    Concepts: matrix vector product^transpose

1164: .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd()
1165: @*/
1166: int MatMultTranspose(Mat mat,Vec x,Vec y)
1167: {
1169:   PetscTruth flg1, flg2;

1174:   MatPreallocated(mat);

1178:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1179:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1180:   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1181: #ifndef PETSC_HAVE_CONSTRAINTS
1182:   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
1183:   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
1184: #endif

1186:   if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP, "Operation not supported");
1187:   PetscLogEventBegin(MAT_MultTranspose,mat,x,y,0);
1188:   if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply tranpose defined");
1189: 
1190:   PetscTypeCompare((PetscObject)mat,MATSEQSBAIJ,&flg1);
1191:   PetscTypeCompare((PetscObject)mat,MATMPISBAIJ,&flg2);
1192:   if (flg1 || flg2) { /* mat is in sbaij format */
1193:     (*mat->ops->mult)(mat,x,y);
1194:   } else {
1195:     (*mat->ops->multtranspose)(mat,x,y);
1196:   }
1197:   PetscLogEventEnd(MAT_MultTranspose,mat,x,y,0);
1198:   return(0);
1199: }

1201: #undef __FUNCT__  
1203: /*@
1204:     MatMultAdd -  Computes v3 = v2 + A * v1.

1206:     Collective on Mat and Vec

1208:     Input Parameters:
1209: +   mat - the matrix
1210: -   v1, v2 - the vectors

1212:     Output Parameters:
1213: .   v3 - the result

1215:     Notes:
1216:     The vectors v1 and v3 cannot be the same.  I.e., one cannot
1217:     call MatMultAdd(A,v1,v2,v1).

1219:     Level: beginner

1221:     Concepts: matrix vector product^addition

1223: .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd()
1224: @*/
1225: int MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1226: {

1232:   MatPreallocated(mat);

1237:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1238:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1239:   if (mat->N != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->N,v1->N);
1240:   if (mat->M != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->M,v2->N);
1241:   if (mat->M != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->M,v3->N);
1242:   if (mat->m != v3->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: local dim %d %d",mat->m,v3->n);
1243:   if (mat->m != v2->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: local dim %d %d",mat->m,v2->n);
1244:   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");

1246:   PetscLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);
1247:   (*mat->ops->multadd)(mat,v1,v2,v3);
1248:   PetscLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);
1249:   return(0);
1250: }

1252: #undef __FUNCT__  
1254: /*@
1255:    MatMultTransposeAdd - Computes v3 = v2 + A' * v1.

1257:    Collective on Mat and Vec

1259:    Input Parameters:
1260: +  mat - the matrix
1261: -  v1, v2 - the vectors

1263:    Output Parameters:
1264: .  v3 - the result

1266:    Notes:
1267:    The vectors v1 and v3 cannot be the same.  I.e., one cannot
1268:    call MatMultTransposeAdd(A,v1,v2,v1).

1270:    Level: beginner

1272:    Concepts: matrix vector product^transpose and addition

1274: .seealso: MatMultTranspose(), MatMultAdd(), MatMult()
1275: @*/
1276: int MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
1277: {

1283:   MatPreallocated(mat);

1288:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1289:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1290:   if (!mat->ops->multtransposeadd) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1291:   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
1292:   if (mat->M != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->M,v1->N);
1293:   if (mat->N != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->N,v2->N);
1294:   if (mat->N != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->N,v3->N);

1296:   PetscLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);
1297:   (*mat->ops->multtransposeadd)(mat,v1,v2,v3);
1298:   PetscLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);
1299:   return(0);
1300: }

1302: #undef __FUNCT__  
1304: /*@
1305:    MatMultConstrained - The inner multiplication routine for a
1306:    constrained matrix P^T A P.

1308:    Collective on Mat and Vec

1310:    Input Parameters:
1311: +  mat - the matrix
1312: -  x   - the vector to be multilplied

1314:    Output Parameters:
1315: .  y - the result

1317:    Notes:
1318:    The vectors x and y cannot be the same.  I.e., one cannot
1319:    call MatMult(A,y,y).

1321:    Level: beginner

1323: .keywords: matrix, multiply, matrix-vector product, constraint
1324: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
1325: @*/
1326: int MatMultConstrained(Mat mat,Vec x,Vec y)
1327: {

1333:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1334:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1335:   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1336:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1337:   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
1338:   if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %d %d",mat->m,y->n);

1340:   PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
1341:   (*mat->ops->multconstrained)(mat,x,y);
1342:   PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);

1344:   return(0);
1345: }

1347: #undef __FUNCT__  
1349: /*@
1350:    MatMultTransposeConstrained - The inner multiplication routine for a
1351:    constrained matrix P^T A^T P.

1353:    Collective on Mat and Vec

1355:    Input Parameters:
1356: +  mat - the matrix
1357: -  x   - the vector to be multilplied

1359:    Output Parameters:
1360: .  y - the result

1362:    Notes:
1363:    The vectors x and y cannot be the same.  I.e., one cannot
1364:    call MatMult(A,y,y).

1366:    Level: beginner

1368: .keywords: matrix, multiply, matrix-vector product, constraint
1369: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
1370: @*/
1371: int MatMultTransposeConstrained(Mat mat,Vec x,Vec y)
1372: {

1378:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1379:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1380:   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
1381:   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1382:   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);

1384:   PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
1385:   (*mat->ops->multtransposeconstrained)(mat,x,y);
1386:   PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);

1388:   return(0);
1389: }
1390: /* ------------------------------------------------------------*/
1391: #undef __FUNCT__  
1393: /*@C
1394:    MatGetInfo - Returns information about matrix storage (number of
1395:    nonzeros, memory, etc.).

1397:    Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used
1398:    as the flag

1400:    Input Parameters:
1401: .  mat - the matrix

1403:    Output Parameters:
1404: +  flag - flag indicating the type of parameters to be returned
1405:    (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
1406:    MAT_GLOBAL_SUM - sum over all processors)
1407: -  info - matrix information context

1409:    Notes:
1410:    The MatInfo context contains a variety of matrix data, including
1411:    number of nonzeros allocated and used, number of mallocs during
1412:    matrix assembly, etc.  Additional information for factored matrices
1413:    is provided (such as the fill ratio, number of mallocs during
1414:    factorization, etc.).  Much of this info is printed to STDOUT
1415:    when using the runtime options 
1416: $       -log_info -mat_view_info

1418:    Example for C/C++ Users:
1419:    See the file ${PETSC_DIR}/include/petscmat.h for a complete list of
1420:    data within the MatInfo context.  For example, 
1421: .vb
1422:       MatInfo info;
1423:       Mat     A;
1424:       double  mal, nz_a, nz_u;

1426:       MatGetInfo(A,MAT_LOCAL,&info);
1427:       mal  = info.mallocs;
1428:       nz_a = info.nz_allocated;
1429: .ve

1431:    Example for Fortran Users:
1432:    Fortran users should declare info as a double precision
1433:    array of dimension MAT_INFO_SIZE, and then extract the parameters
1434:    of interest.  See the file ${PETSC_DIR}/include/finclude/petscmat.h
1435:    a complete list of parameter names.
1436: .vb
1437:       double  precision info(MAT_INFO_SIZE)
1438:       double  precision mal, nz_a
1439:       Mat     A
1440:       integer ierr

1442:       call MatGetInfo(A,MAT_LOCAL,info,ierr)
1443:       mal = info(MAT_INFO_MALLOCS)
1444:       nz_a = info(MAT_INFO_NZ_ALLOCATED)
1445: .ve

1447:     Level: intermediate

1449:     Concepts: matrices^getting information on
1450:  
1451: @*/
1452: int MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
1453: {

1459:   MatPreallocated(mat);
1461:   if (!mat->ops->getinfo) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1462:   (*mat->ops->getinfo)(mat,flag,info);
1463:   return(0);
1464: }

1466: /* ----------------------------------------------------------*/
1467: #undef __FUNCT__  
1469: /*@C  
1470:    MatILUDTFactor - Performs a drop tolerance ILU factorization.

1472:    Collective on Mat

1474:    Input Parameters:
1475: +  mat - the matrix
1476: .  info - information about the factorization to be done
1477: .  row - row permutation
1478: -  col - column permutation

1480:    Output Parameters:
1481: .  fact - the factored matrix

1483:    Level: developer

1485:    Notes:
1486:    Most users should employ the simplified SLES interface for linear solvers
1487:    instead of working directly with matrix algebra routines such as this.
1488:    See, e.g., SLESCreate().

1490:    This is currently only supported for the SeqAIJ matrix format using code
1491:    from Yousef Saad's SPARSEKIT2  package (translated to C with f2c) and/or
1492:    Matlab. SPARSEKIT2 is copyrighted by Yousef Saad with the GNU copyright
1493:    and thus can be distributed with PETSc.

1495:     Concepts: matrices^ILUDT factorization

1497: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatILUInfo
1498: @*/
1499: int MatILUDTFactor(Mat mat,MatILUInfo *info,IS row,IS col,Mat *fact)
1500: {

1506:   MatPreallocated(mat);
1508:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1509:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1510:   if (!mat->ops->iludtfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

1512:   PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
1513:   (*mat->ops->iludtfactor)(mat,info,row,col,fact);
1514:   PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);

1516:   return(0);
1517: }

1519: #undef __FUNCT__  
1521: /*@  
1522:    MatLUFactor - Performs in-place LU factorization of matrix.

1524:    Collective on Mat

1526:    Input Parameters:
1527: +  mat - the matrix
1528: .  row - row permutation
1529: .  col - column permutation
1530: -  info - options for factorization, includes 
1531: $          fill - expected fill as ratio of original fill.
1532: $          dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
1533: $                   Run with the option -log_info to determine an optimal value to use

1535:    Notes:
1536:    Most users should employ the simplified SLES interface for linear solvers
1537:    instead of working directly with matrix algebra routines such as this.
1538:    See, e.g., SLESCreate().

1540:    This changes the state of the matrix to a factored matrix; it cannot be used
1541:    for example with MatSetValues() unless one first calls MatSetUnfactored().

1543:    Level: developer

1545:    Concepts: matrices^LU factorization

1547: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(),
1548:           MatGetOrdering(), MatSetUnfactored(), MatLUInfo

1550: @*/
1551: int MatLUFactor(Mat mat,IS row,IS col,MatLUInfo *info)
1552: {

1558:   MatPreallocated(mat);
1559:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1560:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1561:   if (!mat->ops->lufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

1563:   PetscLogEventBegin(MAT_LUFactor,mat,row,col,0);
1564:   (*mat->ops->lufactor)(mat,row,col,info);
1565:   PetscLogEventEnd(MAT_LUFactor,mat,row,col,0);
1566:   return(0);
1567: }

1569: #undef __FUNCT__  
1571: /*@  
1572:    MatILUFactor - Performs in-place ILU factorization of matrix.

1574:    Collective on Mat

1576:    Input Parameters:
1577: +  mat - the matrix
1578: .  row - row permutation
1579: .  col - column permutation
1580: -  info - structure containing 
1581: $      levels - number of levels of fill.
1582: $      expected fill - as ratio of original fill.
1583: $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
1584:                 missing diagonal entries)

1586:    Notes: 
1587:    Probably really in-place only when level of fill is zero, otherwise allocates
1588:    new space to store factored matrix and deletes previous memory.

1590:    Most users should employ the simplified SLES interface for linear solvers
1591:    instead of working directly with matrix algebra routines such as this.
1592:    See, e.g., SLESCreate().

1594:    Level: developer

1596:    Concepts: matrices^ILU factorization

1598: .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatILUInfo
1599: @*/
1600: int MatILUFactor(Mat mat,IS row,IS col,MatILUInfo *info)
1601: {

1607:   MatPreallocated(mat);
1608:   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
1609:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1610:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1611:   if (!mat->ops->ilufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

1613:   PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
1614:   (*mat->ops->ilufactor)(mat,row,col,info);
1615:   PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);
1616:   return(0);
1617: }

1619: #undef __FUNCT__  
1621: /*@  
1622:    MatLUFactorSymbolic - Performs symbolic LU factorization of matrix.
1623:    Call this routine before calling MatLUFactorNumeric().

1625:    Collective on Mat

1627:    Input Parameters:
1628: +  mat - the matrix
1629: .  row, col - row and column permutations
1630: -  info - options for factorization, includes 
1631: $          fill - expected fill as ratio of original fill.
1632: $          dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
1633: $                   Run with the option -log_info to determine an optimal value to use

1635:    Output Parameter:
1636: .  fact - new matrix that has been symbolically factored

1638:    Notes:
1639:    See the users manual for additional information about
1640:    choosing the fill factor for better efficiency.

1642:    Most users should employ the simplified SLES interface for linear solvers
1643:    instead of working directly with matrix algebra routines such as this.
1644:    See, e.g., SLESCreate().

1646:    Level: developer

1648:    Concepts: matrices^LU symbolic factorization

1650: .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor(), MatLUInfo
1651: @*/
1652: int MatLUFactorSymbolic(Mat mat,IS row,IS col,MatLUInfo *info,Mat *fact)
1653: {

1659:   MatPreallocated(mat);
1663:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1664:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1665:   if (!mat->ops->lufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s  symbolic LU",mat->type_name);

1667:   PetscLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);
1668:   (*mat->ops->lufactorsymbolic)(mat,row,col,info,fact);
1669:   PetscLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);
1670:   return(0);
1671: }

1673: #undef __FUNCT__  
1675: /*@  
1676:    MatLUFactorNumeric - Performs numeric LU factorization of a matrix.
1677:    Call this routine after first calling MatLUFactorSymbolic().

1679:    Collective on Mat

1681:    Input Parameters:
1682: +  mat - the matrix
1683: -  fact - the matrix generated for the factor, from MatLUFactorSymbolic()

1685:    Notes:
1686:    See MatLUFactor() for in-place factorization.  See 
1687:    MatCholeskyFactorNumeric() for the symmetric, positive definite case.  

1689:    Most users should employ the simplified SLES interface for linear solvers
1690:    instead of working directly with matrix algebra routines such as this.
1691:    See, e.g., SLESCreate().

1693:    Level: developer

1695:    Concepts: matrices^LU numeric factorization

1697: .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
1698: @*/
1699: int MatLUFactorNumeric(Mat mat,Mat *fact)
1700: {
1701:   int        ierr;
1702:   PetscTruth flg;

1707:   MatPreallocated(mat);
1710:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1711:   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1712:     SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dimensions are different %d should = %d %d should = %d",
1713:             mat->M,(*fact)->M,mat->N,(*fact)->N);
1714:   }
1715:   if (!(*fact)->ops->lufactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

1717:   PetscLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0);
1718:   (*(*fact)->ops->lufactornumeric)(mat,fact);
1719:   PetscLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0);
1720:   PetscOptionsHasName(PETSC_NULL,"-mat_view_draw",&flg);
1721:   if (flg) {
1722:     PetscOptionsHasName(PETSC_NULL,"-mat_view_contour",&flg);
1723:     if (flg) {
1724:       PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);
1725:     }
1726:     MatView(*fact,PETSC_VIEWER_DRAW_(mat->comm));
1727:     PetscViewerFlush(PETSC_VIEWER_DRAW_(mat->comm));
1728:     if (flg) {
1729:       PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));
1730:     }
1731:   }
1732:   return(0);
1733: }

1735: #undef __FUNCT__  
1737: /*@  
1738:    MatCholeskyFactor - Performs in-place Cholesky factorization of a
1739:    symmetric matrix. 

1741:    Collective on Mat

1743:    Input Parameters:
1744: +  mat - the matrix
1745: .  perm - row and column permutations
1746: -  f - expected fill as ratio of original fill

1748:    Notes:
1749:    See MatLUFactor() for the nonsymmetric case.  See also
1750:    MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().

1752:    Most users should employ the simplified SLES interface for linear solvers
1753:    instead of working directly with matrix algebra routines such as this.
1754:    See, e.g., SLESCreate().

1756:    Level: developer

1758:    Concepts: matrices^Cholesky factorization

1760: .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
1761:           MatGetOrdering()

1763: @*/
1764: int MatCholeskyFactor(Mat mat,IS perm,PetscReal f)
1765: {

1771:   MatPreallocated(mat);
1773:   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
1774:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1775:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1776:   if (!mat->ops->choleskyfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

1778:   PetscLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);
1779:   (*mat->ops->choleskyfactor)(mat,perm,f);
1780:   PetscLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);
1781:   return(0);
1782: }

1784: #undef __FUNCT__  
1786: /*@  
1787:    MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
1788:    of a symmetric matrix. 

1790:    Collective on Mat

1792:    Input Parameters:
1793: +  mat - the matrix
1794: .  perm - row and column permutations
1795: -  f - expected fill as ratio of original

1797:    Output Parameter:
1798: .  fact - the factored matrix

1800:    Notes:
1801:    See MatLUFactorSymbolic() for the nonsymmetric case.  See also
1802:    MatCholeskyFactor() and MatCholeskyFactorNumeric().

1804:    Most users should employ the simplified SLES interface for linear solvers
1805:    instead of working directly with matrix algebra routines such as this.
1806:    See, e.g., SLESCreate().

1808:    Level: developer

1810:    Concepts: matrices^Cholesky symbolic factorization

1812: .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
1813:           MatGetOrdering()

1815: @*/
1816: int MatCholeskyFactorSymbolic(Mat mat,IS perm,PetscReal f,Mat *fact)
1817: {

1823:   MatPreallocated(mat);
1825:   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
1826:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1827:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
1828:   if (!mat->ops->choleskyfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

1830:   PetscLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
1831:   (*mat->ops->choleskyfactorsymbolic)(mat,perm,f,fact);
1832:   PetscLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
1833:   return(0);
1834: }

1836: #undef __FUNCT__  
1838: /*@  
1839:    MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
1840:    of a symmetric matrix. Call this routine after first calling
1841:    MatCholeskyFactorSymbolic().

1843:    Collective on Mat

1845:    Input Parameter:
1846: .  mat - the initial matrix

1848:    Output Parameter:
1849: .  fact - the factored matrix

1851:    Notes:
1852:    Most users should employ the simplified SLES interface for linear solvers
1853:    instead of working directly with matrix algebra routines such as this.
1854:    See, e.g., SLESCreate().

1856:    Level: developer

1858:    Concepts: matrices^Cholesky numeric factorization

1860: .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
1861: @*/
1862: int MatCholeskyFactorNumeric(Mat mat,Mat *fact)
1863: {

1869:   MatPreallocated(mat);
1871:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
1872:   if (!(*fact)->ops->choleskyfactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1873:   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
1874:     SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dim %d should = %d %d should = %d",
1875:             mat->M,(*fact)->M,mat->N,(*fact)->N);
1876:   }

1878:   PetscLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
1879:   (*(*fact)->ops->choleskyfactornumeric)(mat,fact);
1880:   PetscLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
1881:   return(0);
1882: }

1884: /* ----------------------------------------------------------------*/
1885: #undef __FUNCT__  
1887: /*@
1888:    MatSolve - Solves A x = b, given a factored matrix.

1890:    Collective on Mat and Vec

1892:    Input Parameters:
1893: +  mat - the factored matrix
1894: -  b - the right-hand-side vector

1896:    Output Parameter:
1897: .  x - the result vector

1899:    Notes:
1900:    The vectors b and x cannot be the same.  I.e., one cannot
1901:    call MatSolve(A,x,x).

1903:    Notes:
1904:    Most users should employ the simplified SLES interface for linear solvers
1905:    instead of working directly with matrix algebra routines such as this.
1906:    See, e.g., SLESCreate().

1908:    Level: developer

1910:    Concepts: matrices^triangular solves

1912: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd()
1913: @*/
1914: int MatSolve(Mat mat,Vec b,Vec x)
1915: {

1921:   MatPreallocated(mat);
1926:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
1927:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
1928:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1929:   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1930:   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
1931:   if (mat->M == 0 && mat->N == 0) return(0);

1933:   if (!mat->ops->solve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1934:   PetscLogEventBegin(MAT_Solve,mat,b,x,0);
1935:   (*mat->ops->solve)(mat,b,x);
1936:   PetscLogEventEnd(MAT_Solve,mat,b,x,0);
1937:   return(0);
1938: }

1940: #undef __FUNCT__  
1942: /* @
1943:    MatForwardSolve - Solves L x = b, given a factored matrix, A = LU.

1945:    Collective on Mat and Vec

1947:    Input Parameters:
1948: +  mat - the factored matrix
1949: -  b - the right-hand-side vector

1951:    Output Parameter:
1952: .  x - the result vector

1954:    Notes:
1955:    MatSolve() should be used for most applications, as it performs
1956:    a forward solve followed by a backward solve.

1958:    The vectors b and x cannot be the same.  I.e., one cannot
1959:    call MatForwardSolve(A,x,x).

1961:    Most users should employ the simplified SLES interface for linear solvers
1962:    instead of working directly with matrix algebra routines such as this.
1963:    See, e.g., SLESCreate().

1965:    Level: developer

1967:    Concepts: matrices^forward solves

1969: .seealso: MatSolve(), MatBackwardSolve()
1970: @ */
1971: int MatForwardSolve(Mat mat,Vec b,Vec x)
1972: {

1978:   MatPreallocated(mat);
1983:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
1984:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
1985:   if (!mat->ops->forwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
1986:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
1987:   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
1988:   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);

1990:   PetscLogEventBegin(MAT_ForwardSolve,mat,b,x,0);
1991:   (*mat->ops->forwardsolve)(mat,b,x);
1992:   PetscLogEventEnd(MAT_ForwardSolve,mat,b,x,0);
1993:   return(0);
1994: }

1996: #undef __FUNCT__  
1998: /* @
1999:    MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.

2001:    Collective on Mat and Vec

2003:    Input Parameters:
2004: +  mat - the factored matrix
2005: -  b - the right-hand-side vector

2007:    Output Parameter:
2008: .  x - the result vector

2010:    Notes:
2011:    MatSolve() should be used for most applications, as it performs
2012:    a forward solve followed by a backward solve.

2014:    The vectors b and x cannot be the same.  I.e., one cannot
2015:    call MatBackwardSolve(A,x,x).

2017:    Most users should employ the simplified SLES interface for linear solvers
2018:    instead of working directly with matrix algebra routines such as this.
2019:    See, e.g., SLESCreate().

2021:    Level: developer

2023:    Concepts: matrices^backward solves

2025: .seealso: MatSolve(), MatForwardSolve()
2026: @ */
2027: int MatBackwardSolve(Mat mat,Vec b,Vec x)
2028: {

2034:   MatPreallocated(mat);
2039:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2040:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2041:   if (!mat->ops->backwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2042:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
2043:   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
2044:   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);

2046:   PetscLogEventBegin(MAT_BackwardSolve,mat,b,x,0);
2047:   (*mat->ops->backwardsolve)(mat,b,x);
2048:   PetscLogEventEnd(MAT_BackwardSolve,mat,b,x,0);
2049:   return(0);
2050: }

2052: #undef __FUNCT__  
2054: /*@
2055:    MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.

2057:    Collective on Mat and Vec

2059:    Input Parameters:
2060: +  mat - the factored matrix
2061: .  b - the right-hand-side vector
2062: -  y - the vector to be added to 

2064:    Output Parameter:
2065: .  x - the result vector

2067:    Notes:
2068:    The vectors b and x cannot be the same.  I.e., one cannot
2069:    call MatSolveAdd(A,x,y,x).

2071:    Most users should employ the simplified SLES interface for linear solvers
2072:    instead of working directly with matrix algebra routines such as this.
2073:    See, e.g., SLESCreate().

2075:    Level: developer

2077:    Concepts: matrices^triangular solves

2079: .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd()
2080: @*/
2081: int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
2082: {
2083:   PetscScalar one = 1.0;
2084:   Vec    tmp;
2085:   int    ierr;

2090:   MatPreallocated(mat);
2097:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2098:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2099:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
2100:   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
2101:   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N);
2102:   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);
2103:   if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n);

2105:   PetscLogEventBegin(MAT_SolveAdd,mat,b,x,y);
2106:   if (mat->ops->solveadd)  {
2107:     (*mat->ops->solveadd)(mat,b,y,x);
2108:   } else {
2109:     /* do the solve then the add manually */
2110:     if (x != y) {
2111:       MatSolve(mat,b,x);
2112:       VecAXPY(&one,y,x);
2113:     } else {
2114:       VecDuplicate(x,&tmp);
2115:       PetscLogObjectParent(mat,tmp);
2116:       VecCopy(x,tmp);
2117:       MatSolve(mat,b,x);
2118:       VecAXPY(&one,tmp,x);
2119:       VecDestroy(tmp);
2120:     }
2121:   }
2122:   PetscLogEventEnd(MAT_SolveAdd,mat,b,x,y);
2123:   return(0);
2124: }

2126: #undef __FUNCT__  
2128: /*@
2129:    MatSolveTranspose - Solves A' x = b, given a factored matrix.

2131:    Collective on Mat and Vec

2133:    Input Parameters:
2134: +  mat - the factored matrix
2135: -  b - the right-hand-side vector

2137:    Output Parameter:
2138: .  x - the result vector

2140:    Notes:
2141:    The vectors b and x cannot be the same.  I.e., one cannot
2142:    call MatSolveTranspose(A,x,x).

2144:    Most users should employ the simplified SLES interface for linear solvers
2145:    instead of working directly with matrix algebra routines such as this.
2146:    See, e.g., SLESCreate().

2148:    Level: developer

2150:    Concepts: matrices^triangular solves

2152: .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd()
2153: @*/
2154: int MatSolveTranspose(Mat mat,Vec b,Vec x)
2155: {

2161:   MatPreallocated(mat);
2166:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2167:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2168:   if (!mat->ops->solvetranspose) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s",mat->type_name);
2169:   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
2170:   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);

2172:   PetscLogEventBegin(MAT_SolveTranspose,mat,b,x,0);
2173:   (*mat->ops->solvetranspose)(mat,b,x);
2174:   PetscLogEventEnd(MAT_SolveTranspose,mat,b,x,0);
2175:   return(0);
2176: }

2178: #undef __FUNCT__  
2180: /*@
2181:    MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a 
2182:                       factored matrix. 

2184:    Collective on Mat and Vec

2186:    Input Parameters:
2187: +  mat - the factored matrix
2188: .  b - the right-hand-side vector
2189: -  y - the vector to be added to 

2191:    Output Parameter:
2192: .  x - the result vector

2194:    Notes:
2195:    The vectors b and x cannot be the same.  I.e., one cannot
2196:    call MatSolveTransposeAdd(A,x,y,x).

2198:    Most users should employ the simplified SLES interface for linear solvers
2199:    instead of working directly with matrix algebra routines such as this.
2200:    See, e.g., SLESCreate().

2202:    Level: developer

2204:    Concepts: matrices^triangular solves

2206: .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose()
2207: @*/
2208: int MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x)
2209: {
2210:   PetscScalar one = 1.0;
2211:   int         ierr;
2212:   Vec         tmp;

2217:   MatPreallocated(mat);
2224:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
2225:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
2226:   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N);
2227:   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N);
2228:   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->N,y->N);
2229:   if (x->n != y->n)   SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n);

2231:   PetscLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);
2232:   if (mat->ops->solvetransposeadd) {
2233:     (*mat->ops->solvetransposeadd)(mat,b,y,x);
2234:   } else {
2235:     /* do the solve then the add manually */
2236:     if (x != y) {
2237:       MatSolveTranspose(mat,b,x);
2238:       VecAXPY(&one,y,x);
2239:     } else {
2240:       VecDuplicate(x,&tmp);
2241:       PetscLogObjectParent(mat,tmp);
2242:       VecCopy(x,tmp);
2243:       MatSolveTranspose(mat,b,x);
2244:       VecAXPY(&one,tmp,x);
2245:       VecDestroy(tmp);
2246:     }
2247:   }
2248:   PetscLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);
2249:   return(0);
2250: }
2251: /* ----------------------------------------------------------------*/

2253: #undef __FUNCT__  
2255: /*@
2256:    MatRelax - Computes one relaxation sweep.

2258:    Collective on Mat and Vec

2260:    Input Parameters:
2261: +  mat - the matrix
2262: .  b - the right hand side
2263: .  omega - the relaxation factor
2264: .  flag - flag indicating the type of SOR (see below)
2265: .  shift -  diagonal shift
2266: -  its - the number of iterations
2267: -  lits - the number of local iterations 

2269:    Output Parameters:
2270: .  x - the solution (can contain an initial guess)

2272:    SOR Flags:
2273: .     SOR_FORWARD_SWEEP - forward SOR
2274: .     SOR_BACKWARD_SWEEP - backward SOR
2275: .     SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
2276: .     SOR_LOCAL_FORWARD_SWEEP - local forward SOR 
2277: .     SOR_LOCAL_BACKWARD_SWEEP - local forward SOR 
2278: .     SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
2279: .     SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies 
2280:          upper/lower triangular part of matrix to
2281:          vector (with omega)
2282: .     SOR_ZERO_INITIAL_GUESS - zero initial guess

2284:    Notes:
2285:    SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
2286:    SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings
2287:    on each processor. 

2289:    Application programmers will not generally use MatRelax() directly,
2290:    but instead will employ the SLES/PC interface.

2292:    Notes for Advanced Users:
2293:    The flags are implemented as bitwise inclusive or operations.
2294:    For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
2295:    to specify a zero initial guess for SSOR.

2297:    Most users should employ the simplified SLES interface for linear solvers
2298:    instead of working directly with matrix algebra routines such as this.
2299:    See, e.g., SLESCreate().

2301:    Level: developer

2303:    Concepts: matrices^relaxation
2304:    Concepts: matrices^SOR
2305:    Concepts: matrices^Gauss-Seidel

2307: @*/
2308: int MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,int its,int lits,Vec x)
2309: {

2315:   MatPreallocated(mat);
2320:   if (!mat->ops->relax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2321:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2322:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2323:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N);
2324:   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N);
2325:   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n);

2327:   PetscLogEventBegin(MAT_Relax,mat,b,x,0);
2328:   ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,lits,x);
2329:   PetscLogEventEnd(MAT_Relax,mat,b,x,0);
2330:   return(0);
2331: }

2333: #undef __FUNCT__  
2335: /*
2336:       Default matrix copy routine.
2337: */
2338: int MatCopy_Basic(Mat A,Mat B,MatStructure str)
2339: {
2340:   int         ierr,i,rstart,rend,nz,*cwork;
2341:   PetscScalar *vwork;

2344:   MatZeroEntries(B);
2345:   MatGetOwnershipRange(A,&rstart,&rend);
2346:   for (i=rstart; i<rend; i++) {
2347:     MatGetRow(A,i,&nz,&cwork,&vwork);
2348:     MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);
2349:     MatRestoreRow(A,i,&nz,&cwork,&vwork);
2350:   }
2351:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2352:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2353:   return(0);
2354: }

2356: #undef __FUNCT__  
2358: /*@C  
2359:    MatCopy - Copys a matrix to another matrix.

2361:    Collective on Mat

2363:    Input Parameters:
2364: +  A - the matrix
2365: -  str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN

2367:    Output Parameter:
2368: .  B - where the copy is put

2370:    Notes:
2371:    If you use SAME_NONZERO_PATTERN then the two matrices had better have the 
2372:    same nonzero pattern or the routine will crash.

2374:    MatCopy() copies the matrix entries of a matrix to another existing
2375:    matrix (after first zeroing the second matrix).  A related routine is
2376:    MatConvert(), which first creates a new matrix and then copies the data.

2378:    Level: intermediate
2379:    
2380:    Concepts: matrices^copying

2382: .seealso: MatConvert(), MatDuplicate()

2384: @*/
2385: int MatCopy(Mat A,Mat B,MatStructure str)
2386: {

2393:   MatPreallocated(A);
2395:   MatPreallocated(B);
2397:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2398:   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2399:   if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim (%d,%d) (%d,%d)",A->M,B->M,
2400:                                              A->N,B->N);

2402:   PetscLogEventBegin(MAT_Copy,A,B,0,0);
2403:   if (A->ops->copy) {
2404:     (*A->ops->copy)(A,B,str);
2405:   } else { /* generic conversion */
2406:     MatCopy_Basic(A,B,str);
2407:   }
2408:   if (A->mapping) {
2409:     if (B->mapping) {ISLocalToGlobalMappingDestroy(B->mapping);B->mapping = 0;}
2410:     MatSetLocalToGlobalMapping(B,A->mapping);
2411:   }
2412:   if (A->bmapping) {
2413:     if (B->bmapping) {ISLocalToGlobalMappingDestroy(B->bmapping);B->bmapping = 0;}
2414:     MatSetLocalToGlobalMappingBlock(B,A->mapping);
2415:   }
2416:   PetscLogEventEnd(MAT_Copy,A,B,0,0);
2417:   return(0);
2418: }

2420:  #include petscsys.h
2421: PetscTruth MatConvertRegisterAllCalled = PETSC_FALSE;
2422: PetscFList MatConvertList              = 0;

2424: #undef __FUNCT__  
2426: /*@C
2427:     MatConvertRegister - Allows one to register a routine that converts a sparse matrix
2428:         from one format to another.

2430:   Not Collective

2432:   Input Parameters:
2433: +   type - the type of matrix (defined in include/petscmat.h), for example, MATSEQAIJ.
2434: -   Converter - the function that reads the matrix from the binary file.

2436:   Level: developer

2438: .seealso: MatConvertRegisterAll(), MatConvert()

2440: @*/
2441: int MatConvertRegister(char *sname,char *path,char *name,int (*function)(Mat,MatType,Mat*))
2442: {
2443:   int  ierr;
2444:   char fullname[PETSC_MAX_PATH_LEN];

2447:   PetscFListConcat(path,name,fullname);
2448:   PetscFListAdd(&MatConvertList,sname,fullname,(void (*)(void))function);
2449:   return(0);
2450: }

2452: #undef __FUNCT__  
2454: /*@C  
2455:    MatConvert - Converts a matrix to another matrix, either of the same
2456:    or different type.

2458:    Collective on Mat

2460:    Input Parameters:
2461: +  mat - the matrix
2462: -  newtype - new matrix type.  Use MATSAME to create a new matrix of the
2463:    same type as the original matrix.

2465:    Output Parameter:
2466: .  M - pointer to place new matrix

2468:    Notes:
2469:    MatConvert() first creates a new matrix and then copies the data from
2470:    the first matrix.  A related routine is MatCopy(), which copies the matrix
2471:    entries of one matrix to another already existing matrix context.

2473:    Level: intermediate

2475:    Concepts: matrices^converting between storage formats

2477: .seealso: MatCopy(), MatDuplicate()
2478: @*/
2479: int MatConvert(Mat mat,MatType newtype,Mat *M)
2480: {
2481:   int        ierr;
2482:   PetscTruth sametype,issame,flg;
2483:   char       convname[256],mtype[256];

2488:   MatPreallocated(mat);
2490:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2491:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

2493:   PetscOptionsGetString(PETSC_NULL,"-matconvert_type",mtype,256,&flg);
2494:   if (flg) {
2495:     newtype = mtype;
2496:   }
2497:   PetscLogEventBegin(MAT_Convert,mat,0,0,0);
2498: 
2499:   PetscTypeCompare((PetscObject)mat,newtype,&sametype);
2500:   PetscStrcmp(newtype,"same",&issame);
2501:   if ((sametype || issame) && mat->ops->duplicate) {
2502:     (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);
2503:   } else {
2504:     int (*conv)(Mat,MatType,Mat*);
2505:     if (!MatConvertRegisterAllCalled) {
2506:       MatConvertRegisterAll(PETSC_NULL);
2507:     }
2508:     PetscFListFind(mat->comm,MatConvertList,newtype,(void(**)(void))&conv);
2509:     if (conv) {
2510:       (*conv)(mat,newtype,M);
2511:     } else {
2512:       PetscStrcpy(convname,"MatConvert_");
2513:       PetscStrcat(convname,mat->type_name);
2514:       PetscStrcat(convname,"_");
2515:       PetscStrcat(convname,newtype);
2516:       PetscStrcat(convname,"_C");
2517:       PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);
2518:       if (conv) {
2519:         (*conv)(mat,newtype,M);
2520:       } else {
2521:         if (mat->ops->convert) {
2522:           (*mat->ops->convert)(mat,newtype,M);
2523:         } else {
2524:           MatConvert_Basic(mat,newtype,M);
2525:         }
2526:       }
2527:     }
2528:   }
2529:   PetscLogEventEnd(MAT_Convert,mat,0,0,0);
2530:   return(0);
2531: }


2534: #undef __FUNCT__  
2536: /*@C  
2537:    MatDuplicate - Duplicates a matrix including the non-zero structure.

2539:    Collective on Mat

2541:    Input Parameters:
2542: +  mat - the matrix
2543: -  op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero
2544:         values as well or not

2546:    Output Parameter:
2547: .  M - pointer to place new matrix

2549:    Level: intermediate

2551:    Concepts: matrices^duplicating

2553: .seealso: MatCopy(), MatConvert()
2554: @*/
2555: int MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
2556: {

2562:   MatPreallocated(mat);
2564:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2565:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

2567:   *M  = 0;
2568:   PetscLogEventBegin(MAT_Convert,mat,0,0,0);
2569:   if (!mat->ops->duplicate) {
2570:     SETERRQ(PETSC_ERR_SUP,"Not written for this matrix type");
2571:   }
2572:   (*mat->ops->duplicate)(mat,op,M);
2573:   if (mat->mapping) {
2574:     MatSetLocalToGlobalMapping(*M,mat->mapping);
2575:   }
2576:   if (mat->bmapping) {
2577:     MatSetLocalToGlobalMappingBlock(*M,mat->mapping);
2578:   }
2579:   PetscLogEventEnd(MAT_Convert,mat,0,0,0);
2580:   return(0);
2581: }

2583: #undef __FUNCT__  
2585: /*@ 
2586:    MatGetDiagonal - Gets the diagonal of a matrix.

2588:    Collective on Mat and Vec

2590:    Input Parameters:
2591: +  mat - the matrix
2592: -  v - the vector for storing the diagonal

2594:    Output Parameter:
2595: .  v - the diagonal of the matrix

2597:    Notes:
2598:    For the SeqAIJ matrix format, this routine may also be called
2599:    on a LU factored matrix; in that case it routines the reciprocal of 
2600:    the diagonal entries in U. It returns the entries permuted by the 
2601:    row and column permutation used during the symbolic factorization.

2603:    Level: intermediate

2605:    Concepts: matrices^accessing diagonals

2607: .seealso: MatGetRow(), MatGetSubmatrices(), MatGetSubmatrix(), MatGetRowMax()
2608: @*/
2609: int MatGetDiagonal(Mat mat,Vec v)
2610: {

2616:   MatPreallocated(mat);
2619:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2620:   if (!mat->ops->getdiagonal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

2622:   (*mat->ops->getdiagonal)(mat,v);
2623:   return(0);
2624: }

2626: #undef __FUNCT__  
2628: /*@ 
2629:    MatGetRowMax - Gets the maximum value (in absolute value) of each
2630:         row of the matrix

2632:    Collective on Mat and Vec

2634:    Input Parameters:
2635: .  mat - the matrix

2637:    Output Parameter:
2638: .  v - the vector for storing the maximums

2640:    Level: intermediate

2642:    Concepts: matrices^getting row maximums

2644: .seealso: MatGetDiagonal(), MatGetSubmatrices(), MatGetSubmatrix()
2645: @*/
2646: int MatGetRowMax(Mat mat,Vec v)
2647: {

2653:   MatPreallocated(mat);
2656:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2657:   if (!mat->ops->getrowmax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

2659:   (*mat->ops->getrowmax)(mat,v);
2660:   return(0);
2661: }

2663: #undef __FUNCT__  
2665: /*@C
2666:    MatTranspose - Computes an in-place or out-of-place transpose of a matrix.

2668:    Collective on Mat

2670:    Input Parameter:
2671: .  mat - the matrix to transpose

2673:    Output Parameters:
2674: .  B - the transpose (or pass in PETSC_NULL for an in-place transpose)

2676:    Level: intermediate

2678:    Concepts: matrices^transposing

2680: .seealso: MatMultTranspose(), MatMultTransposeAdd()
2681: @*/
2682: int MatTranspose(Mat mat,Mat *B)
2683: {

2689:   MatPreallocated(mat);
2690:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2691:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2692:   if (!mat->ops->transpose) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

2694:   PetscLogEventBegin(MAT_Transpose,mat,0,0,0);
2695:   (*mat->ops->transpose)(mat,B);
2696:   PetscLogEventEnd(MAT_Transpose,mat,0,0,0);
2697:   return(0);
2698: }

2700: #undef __FUNCT__  
2702: /*@C
2703:    MatPermute - Creates a new matrix with rows and columns permuted from the 
2704:    original.

2706:    Collective on Mat

2708:    Input Parameters:
2709: +  mat - the matrix to permute
2710: .  row - row permutation, each processor supplies only the permutation for its rows
2711: -  col - column permutation, each processor needs the entire column permutation, that is
2712:          this is the same size as the total number of columns in the matrix

2714:    Output Parameters:
2715: .  B - the permuted matrix

2717:    Level: advanced

2719:    Concepts: matrices^permuting

2721: .seealso: MatGetOrdering()
2722: @*/
2723: int MatPermute(Mat mat,IS row,IS col,Mat *B)
2724: {

2730:   MatPreallocated(mat);
2733:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2734:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2735:   if (!mat->ops->permute) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2736:   (*mat->ops->permute)(mat,row,col,B);
2737:   return(0);
2738: }

2740: #undef __FUNCT__  
2742: /*@C
2743:   MatPermuteSparsify - Creates a new matrix with rows and columns permuted from the 
2744:   original and sparsified to the prescribed tolerance.

2746:   Collective on Mat

2748:   Input Parameters:
2749: + A    - The matrix to permute
2750: . band - The half-bandwidth of the sparsified matrix, or PETSC_DECIDE
2751: . frac - The half-bandwidth as a fraction of the total size, or 0.0
2752: . tol  - The drop tolerance
2753: . rowp - The row permutation
2754: - colp - The column permutation

2756:   Output Parameter:
2757: . B    - The permuted, sparsified matrix

2759:   Level: advanced

2761:   Note:
2762:   The default behavior (band = PETSC_DECIDE and frac = 0.0) is to
2763:   restrict the half-bandwidth of the resulting matrix to 5% of the
2764:   total matrix size.

2766: .keywords: matrix, permute, sparsify

2768: .seealso: MatGetOrdering(), MatPermute()
2769: @*/
2770: int MatPermuteSparsify(Mat A, int band, PetscReal frac, PetscReal tol, IS rowp, IS colp, Mat *B)
2771: {
2772:   IS           irowp, icolp;
2773:   int         *rows, *cols;
2774:   int          M, N, locRowStart, locRowEnd;
2775:   int          nz, newNz;
2776:   int         *cwork, *cnew;
2777:   PetscScalar *vwork, *vnew;
2778:   int          bw, size;
2779:   int          row, locRow, newRow, col, newCol;
2780:   int          ierr;

2786:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
2787:   if (A->factor)     SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2788:   if (!A->ops->permutesparsify) {
2789:     MatGetSize(A, &M, &N);
2790:     MatGetOwnershipRange(A, &locRowStart, &locRowEnd);
2791:     ISGetSize(rowp, &size);
2792:     if (size != M) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %d for row permutation, should be %d", size, M);
2793:     ISGetSize(colp, &size);
2794:     if (size != N) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %d for column permutation, should be %d", size, N);
2795:     ISInvertPermutation(rowp, 0, &irowp);
2796:     ISGetIndices(irowp, &rows);
2797:     ISInvertPermutation(colp, 0, &icolp);
2798:     ISGetIndices(icolp, &cols);
2799:     PetscMalloc(N * sizeof(int),         &cnew);
2800:     PetscMalloc(N * sizeof(PetscScalar), &vnew);

2802:     /* Setup bandwidth to include */
2803:     if (band == PETSC_DECIDE) {
2804:       if (frac <= 0.0)
2805:         bw = (int) (M * 0.05);
2806:       else
2807:         bw = (int) (M * frac);
2808:     } else {
2809:       if (band <= 0) SETERRQ(PETSC_ERR_ARG_WRONG, "Bandwidth must be a positive integer");
2810:       bw = band;
2811:     }

2813:     /* Put values into new matrix */
2814:     MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, B);
2815:     for(row = locRowStart, locRow = 0; row < locRowEnd; row++, locRow++) {
2816:       MatGetRow(A, row, &nz, &cwork, &vwork);
2817:       newRow   = rows[locRow]+locRowStart;
2818:       for(col = 0, newNz = 0; col < nz; col++) {
2819:         newCol = cols[cwork[col]];
2820:         if ((newCol >= newRow - bw) && (newCol < newRow + bw) && (PetscAbsScalar(vwork[col]) >= tol)) {
2821:           cnew[newNz] = newCol;
2822:           vnew[newNz] = vwork[col];
2823:           newNz++;
2824:         }
2825:       }
2826:       MatSetValues(*B, 1, &newRow, newNz, cnew, vnew, INSERT_VALUES);
2827:       MatRestoreRow(A, row, &nz, &cwork, &vwork);
2828:     }
2829:     PetscFree(cnew);
2830:     PetscFree(vnew);
2831:     MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
2832:     MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
2833:     ISRestoreIndices(irowp, &rows);
2834:     ISRestoreIndices(icolp, &cols);
2835:     ISDestroy(irowp);
2836:     ISDestroy(icolp);
2837:   } else {
2838:     (*A->ops->permutesparsify)(A, band, frac, tol, rowp, colp, B);
2839:   }
2840:   return(0);
2841: }

2843: #undef __FUNCT__  
2845: /*@
2846:    MatEqual - Compares two matrices.

2848:    Collective on Mat

2850:    Input Parameters:
2851: +  A - the first matrix
2852: -  B - the second matrix

2854:    Output Parameter:
2855: .  flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.

2857:    Level: intermediate

2859:    Concepts: matrices^equality between
2860: @*/
2861: int MatEqual(Mat A,Mat B,PetscTruth *flg)
2862: {

2869:   MatPreallocated(A);
2871:   MatPreallocated(B);
2874:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2875:   if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2876:   if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim %d %d %d %d",A->M,B->M,A->N,B->N);
2877:   if (!A->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",A->type_name);
2878:   (*A->ops->equal)(A,B,flg);
2879:   return(0);
2880: }

2882: #undef __FUNCT__  
2884: /*@
2885:    MatDiagonalScale - Scales a matrix on the left and right by diagonal
2886:    matrices that are stored as vectors.  Either of the two scaling
2887:    matrices can be PETSC_NULL.

2889:    Collective on Mat

2891:    Input Parameters:
2892: +  mat - the matrix to be scaled
2893: .  l - the left scaling vector (or PETSC_NULL)
2894: -  r - the right scaling vector (or PETSC_NULL)

2896:    Notes:
2897:    MatDiagonalScale() computes A = LAR, where
2898:    L = a diagonal matrix, R = a diagonal matrix

2900:    Level: intermediate

2902:    Concepts: matrices^diagonal scaling
2903:    Concepts: diagonal scaling of matrices

2905: .seealso: MatScale()
2906: @*/
2907: int MatDiagonalScale(Mat mat,Vec l,Vec r)
2908: {

2914:   MatPreallocated(mat);
2915:   if (!mat->ops->diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2918:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2919:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

2921:   PetscLogEventBegin(MAT_Scale,mat,0,0,0);
2922:   (*mat->ops->diagonalscale)(mat,l,r);
2923:   PetscLogEventEnd(MAT_Scale,mat,0,0,0);
2924:   return(0);
2925: }

2927: #undef __FUNCT__  
2929: /*@
2930:     MatScale - Scales all elements of a matrix by a given number.

2932:     Collective on Mat

2934:     Input Parameters:
2935: +   mat - the matrix to be scaled
2936: -   a  - the scaling value

2938:     Output Parameter:
2939: .   mat - the scaled matrix

2941:     Level: intermediate

2943:     Concepts: matrices^scaling all entries

2945: .seealso: MatDiagonalScale()
2946: @*/
2947: int MatScale(PetscScalar *a,Mat mat)
2948: {

2954:   MatPreallocated(mat);
2956:   if (!mat->ops->scale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2957:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2958:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

2960:   PetscLogEventBegin(MAT_Scale,mat,0,0,0);
2961:   (*mat->ops->scale)(a,mat);
2962:   PetscLogEventEnd(MAT_Scale,mat,0,0,0);
2963:   return(0);
2964: }

2966: #undef __FUNCT__  
2968: /*@ 
2969:    MatNorm - Calculates various norms of a matrix.

2971:    Collective on Mat

2973:    Input Parameters:
2974: +  mat - the matrix
2975: -  type - the type of norm, NORM_1, NORM_FROBENIUS, NORM_INFINITY

2977:    Output Parameters:
2978: .  nrm - the resulting norm 

2980:    Level: intermediate

2982:    Concepts: matrices^norm
2983:    Concepts: norm^of matrix
2984: @*/
2985: int MatNorm(Mat mat,NormType type,PetscReal *nrm)
2986: {

2992:   MatPreallocated(mat);

2995:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
2996:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
2997:   if (!mat->ops->norm) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
2998:   (*mat->ops->norm)(mat,type,nrm);
2999:   return(0);
3000: }

3002: /* 
3003:      This variable is used to prevent counting of MatAssemblyBegin() that
3004:    are called from within a MatAssemblyEnd().
3005: */
3006: static int MatAssemblyEnd_InUse = 0;
3007: #undef __FUNCT__  
3009: /*@
3010:    MatAssemblyBegin - Begins assembling the matrix.  This routine should
3011:    be called after completing all calls to MatSetValues().

3013:    Collective on Mat

3015:    Input Parameters:
3016: +  mat - the matrix 
3017: -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
3018:  
3019:    Notes: 
3020:    MatSetValues() generally caches the values.  The matrix is ready to
3021:    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
3022:    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
3023:    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
3024:    using the matrix.

3026:    Level: beginner

3028:    Concepts: matrices^assembling

3030: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled()
3031: @*/
3032: int MatAssemblyBegin(Mat mat,MatAssemblyType type)
3033: {

3039:   MatPreallocated(mat);
3040:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.nDid you forget to call MatSetUnfactored()?");
3041:   if (mat->assembled) {
3042:     mat->was_assembled = PETSC_TRUE;
3043:     mat->assembled     = PETSC_FALSE;
3044:   }
3045:   if (!MatAssemblyEnd_InUse) {
3046:     PetscLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);
3047:     if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
3048:     PetscLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);
3049:   } else {
3050:     if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
3051:   }
3052:   return(0);
3053: }

3055: #undef __FUNCT__  
3057: /*@
3058:    MatAssembled - Indicates if a matrix has been assembled and is ready for
3059:      use; for example, in matrix-vector product.

3061:    Collective on Mat

3063:    Input Parameter:
3064: .  mat - the matrix 

3066:    Output Parameter:
3067: .  assembled - PETSC_TRUE or PETSC_FALSE

3069:    Level: advanced

3071:    Concepts: matrices^assembled?

3073: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin()
3074: @*/
3075: int MatAssembled(Mat mat,PetscTruth *assembled)
3076: {
3080:   MatPreallocated(mat);
3081:   *assembled = mat->assembled;
3082:   return(0);
3083: }

3085: #undef __FUNCT__  
3087: /*
3088:     Processes command line options to determine if/how a matrix
3089:   is to be viewed. Called by MatAssemblyEnd() and MatLoad().
3090: */
3091: int MatView_Private(Mat mat)
3092: {
3093:   int               ierr;
3094:   PetscTruth        flg;
3095:   static PetscTruth incall = PETSC_FALSE;

3098:   if (incall) return(0);
3099:   incall = PETSC_TRUE;
3100:   PetscOptionsBegin(mat->comm,mat->prefix,"Matrix Options","Mat");
3101:     PetscOptionsName("-mat_view_info","Information on matrix size","MatView",&flg);
3102:     if (flg) {
3103:       PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO);
3104:       MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
3105:       PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));
3106:     }
3107:     PetscOptionsName("-mat_view_info_detailed","Nonzeros in the matrix","MatView",&flg);
3108:     if (flg) {
3109:       PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO_DETAIL);
3110:       MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
3111:       PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));
3112:     }
3113:     PetscOptionsName("-mat_view","Print matrix to stdout","MatView",&flg);
3114:     if (flg) {
3115:       MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
3116:     }
3117:     PetscOptionsName("-mat_view_matlab","Print matrix to stdout in a format Matlab can read","MatView",&flg);
3118:     if (flg) {
3119:       PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_MATLAB);
3120:       MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
3121:       PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));
3122:     }
3123:     PetscOptionsName("-mat_view_socket","Send matrix to socket (can be read from matlab)","MatView",&flg);
3124:     if (flg) {
3125:       MatView(mat,PETSC_VIEWER_SOCKET_(mat->comm));
3126:       PetscViewerFlush(PETSC_VIEWER_SOCKET_(mat->comm));
3127:     }
3128:     PetscOptionsName("-mat_view_binary","Save matrix to file in binary format","MatView",&flg);
3129:     if (flg) {
3130:       MatView(mat,PETSC_VIEWER_BINARY_(mat->comm));
3131:       PetscViewerFlush(PETSC_VIEWER_BINARY_(mat->comm));
3132:     }
3133:   PetscOptionsEnd();
3134:   /* cannot have inside PetscOptionsBegin() because uses PetscOptionsBegin() */
3135:   PetscOptionsHasName(mat->prefix,"-mat_view_draw",&flg);
3136:   if (flg) {
3137:     PetscOptionsHasName(mat->prefix,"-mat_view_contour",&flg);
3138:     if (flg) {
3139:       PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);
3140:     }
3141:     MatView(mat,PETSC_VIEWER_DRAW_(mat->comm));
3142:     PetscViewerFlush(PETSC_VIEWER_DRAW_(mat->comm));
3143:     if (flg) {
3144:       PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));
3145:     }
3146:   }
3147:   incall = PETSC_FALSE;
3148:   return(0);
3149: }

3151: #undef __FUNCT__  
3153: /*@
3154:    MatAssemblyEnd - Completes assembling the matrix.  This routine should
3155:    be called after MatAssemblyBegin().

3157:    Collective on Mat

3159:    Input Parameters:
3160: +  mat - the matrix 
3161: -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY

3163:    Options Database Keys:
3164: +  -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
3165: .  -mat_view_info_detailed - Prints more detailed info
3166: .  -mat_view - Prints matrix in ASCII format
3167: .  -mat_view_matlab - Prints matrix in Matlab format
3168: .  -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
3169: .  -display <name> - Sets display name (default is host)
3170: .  -draw_pause <sec> - Sets number of seconds to pause after display
3171: .  -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
3172: .  -viewer_socket_machine <machine>
3173: .  -viewer_socket_port <port>
3174: .  -mat_view_binary - save matrix to file in binary format
3175: -  -viewer_binary_filename <name>

3177:    Notes: 
3178:    MatSetValues() generally caches the values.  The matrix is ready to
3179:    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
3180:    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
3181:    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
3182:    using the matrix.

3184:    Level: beginner

3186: .seealso: MatAssemblyBegin(), MatSetValues(), PetscDrawOpenX(), MatView(), MatAssembled(), PetscViewerSocketOpen()
3187: @*/
3188: int MatAssemblyEnd(Mat mat,MatAssemblyType type)
3189: {
3190:   int        ierr;
3191:   static int inassm = 0;
3192:   PetscTruth flg;

3197:   MatPreallocated(mat);

3199:   inassm++;
3200:   MatAssemblyEnd_InUse++;
3201:   if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
3202:     PetscLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);
3203:     if (mat->ops->assemblyend) {
3204:       (*mat->ops->assemblyend)(mat,type);
3205:     }
3206:     PetscLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);
3207:   } else {
3208:     if (mat->ops->assemblyend) {
3209:       (*mat->ops->assemblyend)(mat,type);
3210:     }
3211:   }

3213:   /* Flush assembly is not a true assembly */
3214:   if (type != MAT_FLUSH_ASSEMBLY) {
3215:     mat->assembled  = PETSC_TRUE; mat->num_ass++;
3216:   }
3217:   mat->insertmode = NOT_SET_VALUES;
3218:   MatAssemblyEnd_InUse--;

3220:   if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) {
3221:     MatView_Private(mat);
3222:   }
3223:   inassm--;
3224:   PetscOptionsHasName(mat->prefix,"-help",&flg);
3225:   if (flg) {
3226:     MatPrintHelp(mat);
3227:   }
3228:   return(0);
3229: }


3232: #undef __FUNCT__  
3234: /*@
3235:    MatCompress - Tries to store the matrix in as little space as 
3236:    possible.  May fail if memory is already fully used, since it
3237:    tries to allocate new space.

3239:    Collective on Mat

3241:    Input Parameters:
3242: .  mat - the matrix 

3244:    Level: advanced

3246: @*/
3247: int MatCompress(Mat mat)
3248: {

3254:   MatPreallocated(mat);
3255:   if (mat->ops->compress) {(*mat->ops->compress)(mat);}
3256:   return(0);
3257: }

3259: #undef __FUNCT__  
3261: /*@
3262:    MatSetOption - Sets a parameter option for a matrix. Some options
3263:    may be specific to certain storage formats.  Some options
3264:    determine how values will be inserted (or added). Sorted, 
3265:    row-oriented input will generally assemble the fastest. The default
3266:    is row-oriented, nonsorted input. 

3268:    Collective on Mat

3270:    Input Parameters:
3271: +  mat - the matrix 
3272: -  option - the option, one of those listed below (and possibly others),
3273:              e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR

3275:    Options Describing Matrix Structure:
3276: +    MAT_SYMMETRIC - symmetric in terms of both structure and value
3277: -    MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure

3279:    Options For Use with MatSetValues():
3280:    Insert a logically dense subblock, which can be
3281: +    MAT_ROW_ORIENTED - row-oriented (default)
3282: .    MAT_COLUMN_ORIENTED - column-oriented
3283: .    MAT_ROWS_SORTED - sorted by row
3284: .    MAT_ROWS_UNSORTED - not sorted by row (default)
3285: .    MAT_COLUMNS_SORTED - sorted by column
3286: -    MAT_COLUMNS_UNSORTED - not sorted by column (default)

3288:    Not these options reflect the data you pass in with MatSetValues(); it has 
3289:    nothing to do with how the data is stored internally in the matrix 
3290:    data structure.

3292:    When (re)assembling a matrix, we can restrict the input for
3293:    efficiency/debugging purposes.  These options include
3294: +    MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be
3295:         allowed if they generate a new nonzero
3296: .    MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed
3297: .    MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if
3298:          they generate a nonzero in a new diagonal (for block diagonal format only)
3299: .    MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only)
3300: .    MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries
3301: .    MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry
3302: -    MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly

3304:    Notes:
3305:    Some options are relevant only for particular matrix types and
3306:    are thus ignored by others.  Other options are not supported by
3307:    certain matrix types and will generate an error message if set.

3309:    If using a Fortran 77 module to compute a matrix, one may need to 
3310:    use the column-oriented option (or convert to the row-oriented 
3311:    format).  

3313:    MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion 
3314:    that would generate a new entry in the nonzero structure is instead
3315:    ignored.  Thus, if memory has not alredy been allocated for this particular 
3316:    data, then the insertion is ignored. For dense matrices, in which
3317:    the entire array is allocated, no entries are ever ignored. 
3318:    Set after the first MatAssemblyEnd()

3320:    MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion 
3321:    that would generate a new entry in the nonzero structure instead produces 
3322:    an error. (Currently supported for AIJ and BAIJ formats only.)
3323:    This is a useful flag when using SAME_NONZERO_PATTERN in calling
3324:    SLESSetOperators() to ensure that the nonzero pattern truely does 
3325:    remain unchanged. Set after the first MatAssemblyEnd()

3327:    MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion 
3328:    that would generate a new entry that has not been preallocated will 
3329:    instead produce an error. (Currently supported for AIJ and BAIJ formats
3330:    only.) This is a useful flag when debugging matrix memory preallocation.

3332:    MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for 
3333:    other processors should be dropped, rather than stashed.
3334:    This is useful if you know that the "owning" processor is also 
3335:    always generating the correct matrix entries, so that PETSc need
3336:    not transfer duplicate entries generated on another processor.
3337:    
3338:    MAT_USE_HASH_TABLE indicates that a hash table be used to improve the
3339:    searches during matrix assembly. When this flag is set, the hash table
3340:    is created during the first Matrix Assembly. This hash table is
3341:    used the next time through, during MatSetVaules()/MatSetVaulesBlocked()
3342:    to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag 
3343:    should be used with MAT_USE_HASH_TABLE flag. This option is currently
3344:    supported by MATMPIBAIJ format only.

3346:    MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries
3347:    are kept in the nonzero structure

3349:    MAT_IGNORE_ZERO_ENTRIES - when using ADD_VALUES for AIJ matrices this will stop
3350:    zero values from creating a zero location in the matrix

3352:    MAT_USE_INODES - indicates using inode version of the code - works with AIJ and 
3353:    ROWBS matrix types

3355:    MAT_DO_NOT_USE_INODES - indicates not using inode version of the code - works
3356:    with AIJ and ROWBS matrix types

3358:    Level: intermediate

3360:    Concepts: matrices^setting options

3362: @*/
3363: int MatSetOption(Mat mat,MatOption op)
3364: {

3370:   MatPreallocated(mat);
3371:   switch (op) {
3372:   case MAT_SYMMETRIC:
3373:     mat->symmetric              = PETSC_TRUE;
3374:     mat->structurally_symmetric = PETSC_TRUE;
3375:     break;
3376:   case MAT_STRUCTURALLY_SYMMETRIC:
3377:     mat->structurally_symmetric = PETSC_TRUE;
3378:     break;
3379:   default:
3380:     if (mat->ops->setoption) {(*mat->ops->setoption)(mat,op);}
3381:     break;
3382:   }
3383:   return(0);
3384: }

3386: #undef __FUNCT__  
3388: /*@
3389:    MatZeroEntries - Zeros all entries of a matrix.  For sparse matrices
3390:    this routine retains the old nonzero structure.

3392:    Collective on Mat

3394:    Input Parameters:
3395: .  mat - the matrix 

3397:    Level: intermediate

3399:    Concepts: matrices^zeroing

3401: .seealso: MatZeroRows()
3402: @*/
3403: int MatZeroEntries(Mat mat)
3404: {

3410:   MatPreallocated(mat);
3411:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3412:   if (!mat->ops->zeroentries) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

3414:   PetscLogEventBegin(MAT_ZeroEntries,mat,0,0,0);
3415:   (*mat->ops->zeroentries)(mat);
3416:   PetscLogEventEnd(MAT_ZeroEntries,mat,0,0,0);
3417:   return(0);
3418: }

3420: #undef __FUNCT__  
3422: /*@C
3423:    MatZeroRows - Zeros all entries (except possibly the main diagonal)
3424:    of a set of rows of a matrix.

3426:    Collective on Mat

3428:    Input Parameters:
3429: +  mat - the matrix
3430: .  is - index set of rows to remove
3431: -  diag - pointer to value put in all diagonals of eliminated rows.
3432:           Note that diag is not a pointer to an array, but merely a
3433:           pointer to a single value.

3435:    Notes:
3436:    For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
3437:    but does not release memory.  For the dense and block diagonal
3438:    formats this does not alter the nonzero structure.

3440:    If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
3441:    of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
3442:    merely zeroed.

3444:    The user can set a value in the diagonal entry (or for the AIJ and
3445:    row formats can optionally remove the main diagonal entry from the
3446:    nonzero structure as well, by passing a null pointer (PETSC_NULL 
3447:    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).

3449:    For the parallel case, all processes that share the matrix (i.e.,
3450:    those in the communicator used for matrix creation) MUST call this
3451:    routine, regardless of whether any rows being zeroed are owned by
3452:    them.

3454:    For the SBAIJ matrix (since only the upper triangular half of the matrix
3455:    is stored) the effect of this call is to also zero the corresponding 
3456:    column.
3457:   
3458:    Level: intermediate

3460:    Concepts: matrices^zeroing rows

3462: .seealso: MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
3463: @*/
3464: int MatZeroRows(Mat mat,IS is,PetscScalar *diag)
3465: {

3471:   MatPreallocated(mat);
3474:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3475:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3476:   if (!mat->ops->zerorows) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);

3478:   (*mat->ops->zerorows)(mat,is,diag);
3479:   MatView_Private(mat);
3480:   return(0);
3481: }

3483: #undef __FUNCT__  
3485: /*@C 
3486:    MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
3487:    of a set of rows of a matrix; using local numbering of rows.

3489:    Collective on Mat

3491:    Input Parameters:
3492: +  mat - the matrix
3493: .  is - index set of rows to remove
3494: -  diag - pointer to value put in all diagonals of eliminated rows.
3495:           Note that diag is not a pointer to an array, but merely a
3496:           pointer to a single value.

3498:    Notes:
3499:    Before calling MatZeroRowsLocal(), the user must first set the
3500:    local-to-global mapping by calling MatSetLocalToGlobalMapping().

3502:    For the AIJ matrix formats this removes the old nonzero structure,
3503:    but does not release memory.  For the dense and block diagonal
3504:    formats this does not alter the nonzero structure.

3506:    If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure
3507:    of the matrix is not changed (even for AIJ and BAIJ matrices) the values are
3508:    merely zeroed.

3510:    The user can set a value in the diagonal entry (or for the AIJ and
3511:    row formats can optionally remove the main diagonal entry from the
3512:    nonzero structure as well, by passing a null pointer (PETSC_NULL
3513:    in C or PETSC_NULL_SCALAR in Fortran) as the final argument).

3515:    Level: intermediate

3517:    Concepts: matrices^zeroing

3519: .seealso: MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
3520: @*/
3521: int MatZeroRowsLocal(Mat mat,IS is,PetscScalar *diag)
3522: {
3524:   IS  newis;

3529:   MatPreallocated(mat);
3532:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3533:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

3535:   if (mat->ops->zerorowslocal) {
3536:     (*mat->ops->zerorowslocal)(mat,is,diag);
3537:   } else {
3538:     if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Need to provide local to global mapping to matrix first");
3539:     ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);
3540:     (*mat->ops->zerorows)(mat,newis,diag);
3541:     ISDestroy(newis);
3542:   }
3543:   return(0);
3544: }

3546: #undef __FUNCT__  
3548: /*@
3549:    MatGetSize - Returns the numbers of rows and columns in a matrix.

3551:    Not Collective

3553:    Input Parameter:
3554: .  mat - the matrix

3556:    Output Parameters:
3557: +  m - the number of global rows
3558: -  n - the number of global columns

3560:    Level: beginner

3562:    Concepts: matrices^size

3564: .seealso: MatGetLocalSize()
3565: @*/
3566: int MatGetSize(Mat mat,int *m,int* n)
3567: {
3570:   if (m) *m = mat->M;
3571:   if (n) *n = mat->N;
3572:   return(0);
3573: }

3575: #undef __FUNCT__  
3577: /*@
3578:    MatGetLocalSize - Returns the number of rows and columns in a matrix
3579:    stored locally.  This information may be implementation dependent, so
3580:    use with care.

3582:    Not Collective

3584:    Input Parameters:
3585: .  mat - the matrix

3587:    Output Parameters:
3588: +  m - the number of local rows
3589: -  n - the number of local columns

3591:    Level: beginner

3593:    Concepts: matrices^local size

3595: .seealso: MatGetSize()
3596: @*/
3597: int MatGetLocalSize(Mat mat,int *m,int* n)
3598: {
3601:   if (m) *m = mat->m;
3602:   if (n) *n = mat->n;
3603:   return(0);
3604: }

3606: #undef __FUNCT__  
3608: /*@
3609:    MatGetOwnershipRange - Returns the range of matrix rows owned by
3610:    this processor, assuming that the matrix is laid out with the first
3611:    n1 rows on the first processor, the next n2 rows on the second, etc.
3612:    For certain parallel layouts this range may not be well defined.

3614:    Not Collective

3616:    Input Parameters:
3617: .  mat - the matrix

3619:    Output Parameters:
3620: +  m - the global index of the first local row
3621: -  n - one more than the global index of the last local row

3623:    Level: beginner

3625:    Concepts: matrices^row ownership
3626: @*/
3627: int MatGetOwnershipRange(Mat mat,int *m,int* n)
3628: {

3634:   MatPreallocated(mat);
3637:   PetscMapGetLocalRange(mat->rmap,m,n);
3638:   return(0);
3639: }

3641: #undef __FUNCT__  
3643: /*@  
3644:    MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
3645:    Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric() 
3646:    to complete the factorization.

3648:    Collective on Mat

3650:    Input Parameters:
3651: +  mat - the matrix
3652: .  row - row permutation
3653: .  column - column permutation
3654: -  info - structure containing 
3655: $      levels - number of levels of fill.
3656: $      expected fill - as ratio of original fill.
3657: $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
3658:                 missing diagonal entries)

3660:    Output Parameters:
3661: .  fact - new matrix that has been symbolically factored

3663:    Notes:
3664:    See the users manual for additional information about
3665:    choosing the fill factor for better efficiency.

3667:    Most users should employ the simplified SLES interface for linear solvers
3668:    instead of working directly with matrix algebra routines such as this.
3669:    See, e.g., SLESCreate().

3671:    Level: developer

3673:   Concepts: matrices^symbolic LU factorization
3674:   Concepts: matrices^factorization
3675:   Concepts: LU^symbolic factorization

3677: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
3678:           MatGetOrdering(), MatILUInfo

3680: @*/
3681: int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatILUInfo *info,Mat *fact)
3682: {

3688:   MatPreallocated(mat);
3692:   if (info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels of fill negative %d",(int)info->levels);
3693:   if (info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",info->fill);
3694:   if (!mat->ops->ilufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s  symbolic ILU",mat->type_name);
3695:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3696:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

3698:   PetscLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);
3699:   (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);
3700:   PetscLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);
3701:   return(0);
3702: }

3704: #undef __FUNCT__  
3706: /*@  
3707:    MatICCFactorSymbolic - Performs symbolic incomplete
3708:    Cholesky factorization for a symmetric matrix.  Use 
3709:    MatCholeskyFactorNumeric() to complete the factorization.

3711:    Collective on Mat

3713:    Input Parameters:
3714: +  mat - the matrix
3715: .  perm - row and column permutation
3716: .  fill - levels of fill
3717: -  f - expected fill as ratio of original fill

3719:    Output Parameter:
3720: .  fact - the factored matrix

3722:    Notes:
3723:    Currently only no-fill factorization is supported.

3725:    Most users should employ the simplified SLES interface for linear solvers
3726:    instead of working directly with matrix algebra routines such as this.
3727:    See, e.g., SLESCreate().

3729:    Level: developer

3731:   Concepts: matrices^symbolic incomplete Cholesky factorization
3732:   Concepts: matrices^factorization
3733:   Concepts: Cholsky^symbolic factorization

3735: .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor()
3736: @*/
3737: int MatICCFactorSymbolic(Mat mat,IS perm,PetscReal f,int fill,Mat *fact)
3738: {

3744:   MatPreallocated(mat);
3747:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
3748:   if (fill < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Fill negative %d",fill);
3749:   if (f < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",f);
3750:   if (!mat->ops->iccfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s  symbolic ICC",mat->type_name);
3751:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");

3753:   PetscLogEventBegin(MAT_ICCFactorSymbolic,mat,perm,0,0);
3754:   (*mat->ops->iccfactorsymbolic)(mat,perm,f,fill,fact);
3755:   PetscLogEventEnd(MAT_ICCFactorSymbolic,mat,perm,0,0);
3756:   return(0);
3757: }

3759: #undef __FUNCT__  
3761: /*@C
3762:    MatGetArray - Returns a pointer to the element values in the matrix.
3763:    The result of this routine is dependent on the underlying matrix data
3764:    structure, and may not even work for certain matrix types.  You MUST
3765:    call MatRestoreArray() when you no longer need to access the array.

3767:    Not Collective

3769:    Input Parameter:
3770: .  mat - the matrix

3772:    Output Parameter:
3773: .  v - the location of the values


3776:    Fortran Note:
3777:    This routine is used differently from Fortran, e.g.,
3778: .vb
3779:         Mat         mat
3780:         PetscScalar mat_array(1)
3781:         PetscOffset i_mat
3782:         int         ierr
3783:         call MatGetArray(mat,mat_array,i_mat,ierr)

3785:   C  Access first local entry in matrix; note that array is
3786:   C  treated as one dimensional
3787:         value = mat_array(i_mat + 1)

3789:         [... other code ...]
3790:         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3791: .ve

3793:    See the Fortran chapter of the users manual and 
3794:    petsc/src/mat/examples/tests for details.

3796:    Level: advanced

3798:    Concepts: matrices^access array

3800: .seealso: MatRestoreArray(), MatGetArrayF90()
3801: @*/
3802: int MatGetArray(Mat mat,PetscScalar **v)
3803: {

3809:   MatPreallocated(mat);
3811:   if (!mat->ops->getarray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3812:   (*mat->ops->getarray)(mat,v);
3813:   return(0);
3814: }

3816: #undef __FUNCT__  
3818: /*@C
3819:    MatRestoreArray - Restores the matrix after MatGetArray() has been called.

3821:    Not Collective

3823:    Input Parameter:
3824: +  mat - the matrix
3825: -  v - the location of the values

3827:    Fortran Note:
3828:    This routine is used differently from Fortran, e.g.,
3829: .vb
3830:         Mat         mat
3831:         PetscScalar mat_array(1)
3832:         PetscOffset i_mat
3833:         int         ierr
3834:         call MatGetArray(mat,mat_array,i_mat,ierr)

3836:   C  Access first local entry in matrix; note that array is
3837:   C  treated as one dimensional
3838:         value = mat_array(i_mat + 1)

3840:         [... other code ...]
3841:         call MatRestoreArray(mat,mat_array,i_mat,ierr)
3842: .ve

3844:    See the Fortran chapter of the users manual and 
3845:    petsc/src/mat/examples/tests for details

3847:    Level: advanced

3849: .seealso: MatGetArray(), MatRestoreArrayF90()
3850: @*/
3851: int MatRestoreArray(Mat mat,PetscScalar **v)
3852: {

3858:   MatPreallocated(mat);
3860:   if (!mat->ops->restorearray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3861:   (*mat->ops->restorearray)(mat,v);
3862:   return(0);
3863: }

3865: #undef __FUNCT__  
3867: /*@C
3868:    MatGetSubMatrices - Extracts several submatrices from a matrix. If submat
3869:    points to an array of valid matrices, they may be reused to store the new
3870:    submatrices.

3872:    Collective on Mat

3874:    Input Parameters:
3875: +  mat - the matrix
3876: .  n   - the number of submatrixes to be extracted (on this processor, may be zero)
3877: .  irow, icol - index sets of rows and columns to extract
3878: -  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

3880:    Output Parameter:
3881: .  submat - the array of submatrices

3883:    Notes:
3884:    MatGetSubMatrices() can extract only sequential submatrices
3885:    (from both sequential and parallel matrices). Use MatGetSubMatrix()
3886:    to extract a parallel submatrix.

3888:    When extracting submatrices from a parallel matrix, each processor can
3889:    form a different submatrix by setting the rows and columns of its
3890:    individual index sets according to the local submatrix desired.

3892:    When finished using the submatrices, the user should destroy
3893:    them with MatDestroyMatrices().

3895:    MAT_REUSE_MATRIX can only be used when the nonzero structure of the 
3896:    original matrix has not changed from that last call to MatGetSubMatrices().

3898:    This routine creates the matrices in submat; you should NOT create them before
3899:    calling it. It also allocates the array of matrix pointers submat.

3901:    Fortran Note:
3902:    The Fortran interface is slightly different from that given below; it 
3903:    requires one to pass in  as submat a Mat (integer) array of size at least m.

3905:    Level: advanced

3907:    Concepts: matrices^accessing submatrices
3908:    Concepts: submatrices

3910: .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal()
3911: @*/
3912: int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatReuse scall,Mat **submat)
3913: {
3914:   int        ierr;

3919:   MatPreallocated(mat);
3920:   if (!mat->ops->getsubmatrices) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3921:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");

3923:   PetscLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);
3924:   (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);
3925:   PetscLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);
3926:   return(0);
3927: }

3929: #undef __FUNCT__  
3931: /*@C
3932:    MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices().

3934:    Collective on Mat

3936:    Input Parameters:
3937: +  n - the number of local matrices
3938: -  mat - the matrices

3940:    Level: advanced

3942:     Notes: Frees not only the matrices, but also the array that contains the matrices

3944: .seealso: MatGetSubMatrices()
3945: @*/
3946: int MatDestroyMatrices(int n,Mat **mat)
3947: {
3948:   int ierr,i;

3951:   if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of matrices %d",n);
3953:   for (i=0; i<n; i++) {
3954:     MatDestroy((*mat)[i]);
3955:   }
3956:   /* memory is allocated even if n = 0 */
3957:   PetscFree(*mat);
3958:   return(0);
3959: }

3961: #undef __FUNCT__  
3963: /*@
3964:    MatIncreaseOverlap - Given a set of submatrices indicated by index sets,
3965:    replaces the index sets by larger ones that represent submatrices with
3966:    additional overlap.

3968:    Collective on Mat

3970:    Input Parameters:
3971: +  mat - the matrix
3972: .  n   - the number of index sets
3973: .  is  - the array of pointers to index sets
3974: -  ov  - the additional overlap requested

3976:    Level: developer

3978:    Concepts: overlap
3979:    Concepts: ASM^computing overlap

3981: .seealso: MatGetSubMatrices()
3982: @*/
3983: int MatIncreaseOverlap(Mat mat,int n,IS *is,int ov)
3984: {

3990:   MatPreallocated(mat);
3991:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
3992:   if (mat->factor)     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

3994:   if (!ov) return(0);
3995:   if (!mat->ops->increaseoverlap) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
3996:   PetscLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);
3997:   (*mat->ops->increaseoverlap)(mat,n,is,ov);
3998:   PetscLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);
3999:   return(0);
4000: }

4002: #undef __FUNCT__  
4004: /*@
4005:    MatPrintHelp - Prints all the options for the matrix.

4007:    Collective on Mat

4009:    Input Parameter:
4010: .  mat - the matrix 

4012:    Options Database Keys:
4013: +  -help - Prints matrix options
4014: -  -h - Prints matrix options

4016:    Level: developer

4018: .seealso: MatCreate(), MatCreateXXX()
4019: @*/
4020: int MatPrintHelp(Mat mat)
4021: {
4022:   static PetscTruth called = PETSC_FALSE;
4023:   int               ierr;

4028:   MatPreallocated(mat);

4030:   if (!called) {
4031:     if (mat->ops->printhelp) {
4032:       (*mat->ops->printhelp)(mat);
4033:     }
4034:     called = PETSC_TRUE;
4035:   }
4036:   return(0);
4037: }

4039: #undef __FUNCT__  
4041: /*@
4042:    MatGetBlockSize - Returns the matrix block size; useful especially for the
4043:    block row and block diagonal formats.
4044:    
4045:    Not Collective

4047:    Input Parameter:
4048: .  mat - the matrix

4050:    Output Parameter:
4051: .  bs - block size

4053:    Notes:
4054:    Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG.
4055:    Block row formats are MATSEQBAIJ, MATMPIBAIJ

4057:    Level: intermediate

4059:    Concepts: matrices^block size

4061: .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag()
4062: @*/
4063: int MatGetBlockSize(Mat mat,int *bs)
4064: {

4070:   MatPreallocated(mat);
4072:   if (!mat->ops->getblocksize) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4073:   (*mat->ops->getblocksize)(mat,bs);
4074:   return(0);
4075: }

4077: #undef __FUNCT__  
4079: /*@C
4080:     MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices.

4082:    Collective on Mat

4084:     Input Parameters:
4085: +   mat - the matrix
4086: .   shift -  0 or 1 indicating we want the indices starting at 0 or 1
4087: -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
4088:                 symmetrized

4090:     Output Parameters:
4091: +   n - number of rows in the (possibly compressed) matrix
4092: .   ia - the row pointers
4093: .   ja - the column indices
4094: -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned

4096:     Level: developer

4098: .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
4099: @*/
4100: int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
4101: {

4107:   MatPreallocated(mat);
4111:   if (!mat->ops->getrowij) *done = PETSC_FALSE;
4112:   else {
4113:     *done = PETSC_TRUE;
4114:     ierr  = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);
4115:   }
4116:   return(0);
4117: }

4119: #undef __FUNCT__  
4121: /*@C
4122:     MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices.

4124:     Collective on Mat

4126:     Input Parameters:
4127: +   mat - the matrix
4128: .   shift - 1 or zero indicating we want the indices starting at 0 or 1
4129: -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
4130:                 symmetrized

4132:     Output Parameters:
4133: +   n - number of columns in the (possibly compressed) matrix
4134: .   ia - the column pointers
4135: .   ja - the row indices
4136: -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned

4138:     Level: developer

4140: .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
4141: @*/
4142: int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
4143: {

4149:   MatPreallocated(mat);

4154:   if (!mat->ops->getcolumnij) *done = PETSC_FALSE;
4155:   else {
4156:     *done = PETSC_TRUE;
4157:     ierr  = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);
4158:   }
4159:   return(0);
4160: }

4162: #undef __FUNCT__  
4164: /*@C
4165:     MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with
4166:     MatGetRowIJ().

4168:     Collective on Mat

4170:     Input Parameters:
4171: +   mat - the matrix
4172: .   shift - 1 or zero indicating we want the indices starting at 0 or 1
4173: -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
4174:                 symmetrized

4176:     Output Parameters:
4177: +   n - size of (possibly compressed) matrix
4178: .   ia - the row pointers
4179: .   ja - the column indices
4180: -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned

4182:     Level: developer

4184: .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
4185: @*/
4186: int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
4187: {

4193:   MatPreallocated(mat);

4198:   if (!mat->ops->restorerowij) *done = PETSC_FALSE;
4199:   else {
4200:     *done = PETSC_TRUE;
4201:     ierr  = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);
4202:   }
4203:   return(0);
4204: }

4206: #undef __FUNCT__  
4208: /*@C
4209:     MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with
4210:     MatGetColumnIJ().

4212:     Collective on Mat

4214:     Input Parameters:
4215: +   mat - the matrix
4216: .   shift - 1 or zero indicating we want the indices starting at 0 or 1
4217: -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
4218:                 symmetrized

4220:     Output Parameters:
4221: +   n - size of (possibly compressed) matrix
4222: .   ia - the column pointers
4223: .   ja - the row indices
4224: -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned

4226:     Level: developer

4228: .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
4229: @*/
4230: int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done)
4231: {

4237:   MatPreallocated(mat);

4242:   if (!mat->ops->restorecolumnij) *done = PETSC_FALSE;
4243:   else {
4244:     *done = PETSC_TRUE;
4245:     ierr  = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);
4246:   }
4247:   return(0);
4248: }

4250: #undef __FUNCT__  
4252: /*@C
4253:     MatColoringPatch -Used inside matrix coloring routines that 
4254:     use MatGetRowIJ() and/or MatGetColumnIJ().

4256:     Collective on Mat

4258:     Input Parameters:
4259: +   mat - the matrix
4260: .   n   - number of colors
4261: -   colorarray - array indicating color for each column

4263:     Output Parameters:
4264: .   iscoloring - coloring generated using colorarray information

4266:     Level: developer

4268: .seealso: MatGetRowIJ(), MatGetColumnIJ()

4270: @*/
4271: int MatColoringPatch(Mat mat,int n,int ncolors,int *colorarray,ISColoring *iscoloring)
4272: {

4278:   MatPreallocated(mat);

4281:   if (!mat->ops->coloringpatch){
4282:     ISColoringCreate(mat->comm,n,colorarray,iscoloring);
4283:   } else {
4284:     (*mat->ops->coloringpatch)(mat,n,ncolors,colorarray,iscoloring);
4285:   }
4286:   return(0);
4287: }


4290: #undef __FUNCT__  
4292: /*@
4293:    MatSetUnfactored - Resets a factored matrix to be treated as unfactored.

4295:    Collective on Mat

4297:    Input Parameter:
4298: .  mat - the factored matrix to be reset

4300:    Notes: 
4301:    This routine should be used only with factored matrices formed by in-place
4302:    factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE
4303:    format).  This option can save memory, for example, when solving nonlinear
4304:    systems with a matrix-free Newton-Krylov method and a matrix-based, in-place
4305:    ILU(0) preconditioner.  

4307:    Note that one can specify in-place ILU(0) factorization by calling 
4308: .vb
4309:      PCType(pc,PCILU);
4310:      PCILUSeUseInPlace(pc);
4311: .ve
4312:    or by using the options -pc_type ilu -pc_ilu_in_place

4314:    In-place factorization ILU(0) can also be used as a local
4315:    solver for the blocks within the block Jacobi or additive Schwarz
4316:    methods (runtime option: -sub_pc_ilu_in_place).  See the discussion 
4317:    of these preconditioners in the users manual for details on setting
4318:    local solver options.

4320:    Most users should employ the simplified SLES interface for linear solvers
4321:    instead of working directly with matrix algebra routines such as this.
4322:    See, e.g., SLESCreate().

4324:    Level: developer

4326: .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace()

4328:    Concepts: matrices^unfactored

4330: @*/
4331: int MatSetUnfactored(Mat mat)
4332: {

4338:   MatPreallocated(mat);
4339:   mat->factor = 0;
4340:   if (!mat->ops->setunfactored) return(0);
4341:   (*mat->ops->setunfactored)(mat);
4342:   return(0);
4343: }

4345: /*MC
4346:     MatGetArrayF90 - Accesses a matrix array from Fortran90.

4348:     Synopsis:
4349:     MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)

4351:     Not collective

4353:     Input Parameter:
4354: .   x - matrix

4356:     Output Parameters:
4357: +   xx_v - the Fortran90 pointer to the array
4358: -   ierr - error code

4360:     Example of Usage: 
4361: .vb
4362:       PetscScalar, pointer xx_v(:)
4363:       ....
4364:       call MatGetArrayF90(x,xx_v,ierr)
4365:       a = xx_v(3)
4366:       call MatRestoreArrayF90(x,xx_v,ierr)
4367: .ve

4369:     Notes:
4370:     Not yet supported for all F90 compilers

4372:     Level: advanced

4374: .seealso:  MatRestoreArrayF90(), MatGetArray(), MatRestoreArray()

4376:     Concepts: matrices^accessing array

4378: M*/

4380: /*MC
4381:     MatRestoreArrayF90 - Restores a matrix array that has been
4382:     accessed with MatGetArrayF90().

4384:     Synopsis:
4385:     MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)

4387:     Not collective

4389:     Input Parameters:
4390: +   x - matrix
4391: -   xx_v - the Fortran90 pointer to the array

4393:     Output Parameter:
4394: .   ierr - error code

4396:     Example of Usage: 
4397: .vb
4398:        PetscScalar, pointer xx_v(:)
4399:        ....
4400:        call MatGetArrayF90(x,xx_v,ierr)
4401:        a = xx_v(3)
4402:        call MatRestoreArrayF90(x,xx_v,ierr)
4403: .ve
4404:    
4405:     Notes:
4406:     Not yet supported for all F90 compilers

4408:     Level: advanced

4410: .seealso:  MatGetArrayF90(), MatGetArray(), MatRestoreArray()

4412: M*/


4415: #undef __FUNCT__  
4417: /*@
4418:     MatGetSubMatrix - Gets a single submatrix on the same number of processors
4419:                       as the original matrix.

4421:     Collective on Mat

4423:     Input Parameters:
4424: +   mat - the original matrix
4425: .   isrow - rows this processor should obtain
4426: .   iscol - columns for all processors you wish to keep
4427: .   csize - number of columns "local" to this processor (does nothing for sequential 
4428:             matrices). This should match the result from VecGetLocalSize(x,...) if you 
4429:             plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE
4430: -   cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

4432:     Output Parameter:
4433: .   newmat - the new submatrix, of the same type as the old

4435:     Level: advanced

4437:     Notes: the iscol argument MUST be the same on each processor. You might be 
4438:     able to create the iscol argument with ISAllGather().

4440:       The first time this is called you should use a cll of MAT_INITIAL_MATRIX,
4441:    the MatGetSubMatrix() routine will create the newmat for you. Any additional calls
4442:    to this routine with a mat of the same nonzero structure will reuse the matrix
4443:    generated the first time.

4445:     Concepts: matrices^submatrices

4447: .seealso: MatGetSubMatrices(), ISAllGather()
4448: @*/
4449: int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat)
4450: {
4451:   int     ierr, size;
4452:   Mat     *local;

4456:   MatPreallocated(mat);
4457:   MPI_Comm_size(mat->comm,&size);

4459:   /* if original matrix is on just one processor then use submatrix generated */
4460:   if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) {
4461:     MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);
4462:     return(0);
4463:   } else if (!mat->ops->getsubmatrix && size == 1) {
4464:     ierr    = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);
4465:     *newmat = *local;
4466:     ierr    = PetscFree(local);
4467:     return(0);
4468:   }

4470:   if (!mat->ops->getsubmatrix) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4471:   (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);
4472:   return(0);
4473: }

4475: #undef __FUNCT__  
4477: /*@C
4478:    MatGetPetscMaps - Returns the maps associated with the matrix.

4480:    Not Collective

4482:    Input Parameter:
4483: .  mat - the matrix

4485:    Output Parameters:
4486: +  rmap - the row (right) map
4487: -  cmap - the column (left) map  

4489:    Level: developer

4491:    Concepts: maps^getting from matrix

4493: @*/
4494: int MatGetPetscMaps(Mat mat,PetscMap *rmap,PetscMap *cmap)
4495: {

4501:   MatPreallocated(mat);
4502:   (*mat->ops->getmaps)(mat,rmap,cmap);
4503:   return(0);
4504: }

4506: /*
4507:       Version that works for all PETSc matrices
4508: */
4509: #undef __FUNCT__  
4511: int MatGetPetscMaps_Petsc(Mat mat,PetscMap *rmap,PetscMap *cmap)
4512: {
4514:   if (rmap) *rmap = mat->rmap;
4515:   if (cmap) *cmap = mat->cmap;
4516:   return(0);
4517: }

4519: #undef __FUNCT__  
4521: /*@
4522:    MatSetStashInitialSize - sets the sizes of the matrix stash, that is
4523:    used during the assembly process to store values that belong to 
4524:    other processors.

4526:    Not Collective

4528:    Input Parameters:
4529: +  mat   - the matrix
4530: .  size  - the initial size of the stash.
4531: -  bsize - the initial size of the block-stash(if used).

4533:    Options Database Keys:
4534: +   -matstash_initial_size <size> or <size0,size1,...sizep-1>
4535: -   -matstash_block_initial_size <bsize>  or <bsize0,bsize1,...bsizep-1>

4537:    Level: intermediate

4539:    Notes: 
4540:      The block-stash is used for values set with VecSetValuesBlocked() while
4541:      the stash is used for values set with VecSetValues()

4543:      Run with the option -log_info and look for output of the form
4544:      MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
4545:      to determine the appropriate value, MM, to use for size and 
4546:      MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
4547:      to determine the value, BMM to use for bsize

4549:    Concepts: stash^setting matrix size
4550:    Concepts: matrices^stash

4552: @*/
4553: int MatSetStashInitialSize(Mat mat,int size, int bsize)
4554: {

4560:   MatPreallocated(mat);
4561:   MatStashSetInitialSize_Private(&mat->stash,size);
4562:   MatStashSetInitialSize_Private(&mat->bstash,bsize);
4563:   return(0);
4564: }

4566: #undef __FUNCT__  
4568: /*@
4569:    MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of 
4570:      the matrix

4572:    Collective on Mat

4574:    Input Parameters:
4575: +  mat   - the matrix
4576: .  x,y - the vectors
4577: -  w - where the result is stored

4579:    Level: intermediate

4581:    Notes: 
4582:     w may be the same vector as y. 

4584:     This allows one to use either the restriction or interpolation (its transpose)
4585:     matrix to do the interpolation

4587:     Concepts: interpolation

4589: .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()

4591: @*/
4592: int MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w)
4593: {
4594:   int M,N,ierr;

4598:   MatPreallocated(A);
4599:   MatGetSize(A,&M,&N);
4600:   if (N > M) {
4601:     MatMultTransposeAdd(A,x,y,w);
4602:   } else {
4603:     MatMultAdd(A,x,y,w);
4604:   }
4605:   return(0);
4606: }

4608: #undef __FUNCT__  
4610: /*@
4611:    MatInterpolate - y = A*x or A'*x depending on the shape of 
4612:      the matrix

4614:    Collective on Mat

4616:    Input Parameters:
4617: +  mat   - the matrix
4618: -  x,y - the vectors

4620:    Level: intermediate

4622:    Notes: 
4623:     This allows one to use either the restriction or interpolation (its transpose)
4624:     matrix to do the interpolation

4626:    Concepts: matrices^interpolation

4628: .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict()

4630: @*/
4631: int MatInterpolate(Mat A,Vec x,Vec y)
4632: {
4633:   int M,N,ierr;

4637:   MatPreallocated(A);
4638:   MatGetSize(A,&M,&N);
4639:   if (N > M) {
4640:     MatMultTranspose(A,x,y);
4641:   } else {
4642:     MatMult(A,x,y);
4643:   }
4644:   return(0);
4645: }

4647: #undef __FUNCT__  
4649: /*@
4650:    MatRestrict - y = A*x or A'*x

4652:    Collective on Mat

4654:    Input Parameters:
4655: +  mat   - the matrix
4656: -  x,y - the vectors

4658:    Level: intermediate

4660:    Notes: 
4661:     This allows one to use either the restriction or interpolation (its transpose)
4662:     matrix to do the restriction

4664:    Concepts: matrices^restriction

4666: .seealso: MatMultAdd(), MatMultTransposeAdd(), MatInterpolate()

4668: @*/
4669: int MatRestrict(Mat A,Vec x,Vec y)
4670: {
4671:   int M,N,ierr;

4675:   MatPreallocated(A);
4676:   MatGetSize(A,&M,&N);
4677:   if (N > M) {
4678:     MatMult(A,x,y);
4679:   } else {
4680:     MatMultTranspose(A,x,y);
4681:   }
4682:   return(0);
4683: }

4685: #undef __FUNCT__  
4687: /*@C
4688:    MatNullSpaceAttach - attaches a null space to a matrix.
4689:         This null space will be removed from the resulting vector whenever
4690:         MatMult() is called

4692:    Collective on Mat

4694:    Input Parameters:
4695: +  mat - the matrix
4696: -  nullsp - the null space object

4698:    Level: developer

4700:    Notes:
4701:       Overwrites any previous null space that may have been attached

4703:    Concepts: null space^attaching to matrix

4705: .seealso: MatCreate(), MatNullSpaceCreate()
4706: @*/
4707: int MatNullSpaceAttach(Mat mat,MatNullSpace nullsp)
4708: {

4714:   MatPreallocated(mat);

4717:   if (mat->nullsp) {
4718:     MatNullSpaceDestroy(mat->nullsp);
4719:   }
4720:   mat->nullsp = nullsp;
4721:   PetscObjectReference((PetscObject)nullsp);
4722:   return(0);
4723: }

4725: #undef __FUNCT__  
4727: /*@  
4728:    MatICCFactor - Performs in-place incomplete Cholesky factorization of matrix.

4730:    Collective on Mat

4732:    Input Parameters:
4733: +  mat - the matrix
4734: .  row - row/column permutation
4735: .  fill - expected fill factor >= 1.0
4736: -  level - level of fill, for ICC(k)

4738:    Notes: 
4739:    Probably really in-place only when level of fill is zero, otherwise allocates
4740:    new space to store factored matrix and deletes previous memory.

4742:    Most users should employ the simplified SLES interface for linear solvers
4743:    instead of working directly with matrix algebra routines such as this.
4744:    See, e.g., SLESCreate().

4746:    Level: developer

4748:    Concepts: matrices^incomplete Cholesky factorization
4749:    Concepts: Cholesky factorization

4751: .seealso: MatICCFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
4752: @*/
4753: int MatICCFactor(Mat mat,IS row,PetscReal fill,int level)
4754: {

4760:   MatPreallocated(mat);
4761:   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
4762:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
4763:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
4764:   if (!mat->ops->iccfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4765:   (*mat->ops->iccfactor)(mat,row,fill,level);
4766:   return(0);
4767: }

4769: #undef __FUNCT__  
4771: /*@ 
4772:    MatSetValuesAdic - Sets values computed with ADIC automatic differentiation into a matrix.

4774:    Not Collective

4776:    Input Parameters:
4777: +  mat - the matrix
4778: -  v - the values compute with ADIC

4780:    Level: developer

4782:    Notes:
4783:      Must call MatSetColoring() before using this routine. Also this matrix must already
4784:      have its nonzero pattern determined.

4786: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
4787:           MatSetValues(), MatSetColoring(), MatSetValuesAdifor()
4788: @*/
4789: int MatSetValuesAdic(Mat mat,void *v)
4790: {


4797:   if (!mat->assembled) {
4798:     SETERRQ(1,"Matrix must be already assembled");
4799:   }
4800:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
4801:   if (!mat->ops->setvaluesadic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4802:   (*mat->ops->setvaluesadic)(mat,v);
4803:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
4804:   MatView_Private(mat);
4805:   return(0);
4806: }


4809: #undef __FUNCT__  
4811: /*@ 
4812:    MatSetColoring - Sets a coloring used by calls to MatSetValuesAdic()

4814:    Not Collective

4816:    Input Parameters:
4817: +  mat - the matrix
4818: -  coloring - the coloring

4820:    Level: developer

4822: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
4823:           MatSetValues(), MatSetValuesAdic()
4824: @*/
4825: int MatSetColoring(Mat mat,ISColoring coloring)
4826: {


4833:   if (!mat->assembled) {
4834:     SETERRQ(1,"Matrix must be already assembled");
4835:   }
4836:   if (!mat->ops->setcoloring) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4837:   (*mat->ops->setcoloring)(mat,coloring);
4838:   return(0);
4839: }

4841: #undef __FUNCT__  
4843: /*@ 
4844:    MatSetValuesAdifor - Sets values computed with automatic differentiation into a matrix.

4846:    Not Collective

4848:    Input Parameters:
4849: +  mat - the matrix
4850: .  nl - leading dimension of v
4851: -  v - the values compute with ADIFOR

4853:    Level: developer

4855:    Notes:
4856:      Must call MatSetColoring() before using this routine. Also this matrix must already
4857:      have its nonzero pattern determined.

4859: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
4860:           MatSetValues(), MatSetColoring()
4861: @*/
4862: int MatSetValuesAdifor(Mat mat,int nl,void *v)
4863: {


4870:   if (!mat->assembled) {
4871:     SETERRQ(1,"Matrix must be already assembled");
4872:   }
4873:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
4874:   if (!mat->ops->setvaluesadifor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4875:   (*mat->ops->setvaluesadifor)(mat,nl,v);
4876:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
4877:   return(0);
4878: }

4880: EXTERN int MatMPIAIJDiagonalScaleLocal(Mat,Vec);
4881: EXTERN int MatMPIBAIJDiagonalScaleLocal(Mat,Vec);

4883: #undef __FUNCT__  
4885: /*@ 
4886:    MatDiagonalScaleLocal - Scales columns of a matrix given the scaling values including the 
4887:          ghosted ones.

4889:    Not Collective

4891:    Input Parameters:
4892: +  mat - the matrix
4893: -  diag = the diagonal values, including ghost ones

4895:    Level: developer

4897:    Notes: Works only for MPIAIJ and MPIBAIJ matrices
4898:       
4899: .seealso: MatDiagonalScale()
4900: @*/
4901: int MatDiagonalScaleLocal(Mat mat,Vec diag)
4902: {
4903:   int        ierr;
4904:   PetscTruth flag;


4911:   if (!mat->assembled) {
4912:     SETERRQ(1,"Matrix must be already assembled");
4913:   }
4914:   PetscLogEventBegin(MAT_Scale,mat,0,0,0);
4915:   PetscTypeCompare((PetscObject)mat,MATMPIAIJ,&flag);
4916:   if (flag) {
4917:     MatMPIAIJDiagonalScaleLocal(mat,diag);
4918:   } else {
4919:     PetscTypeCompare((PetscObject)mat,MATMPIBAIJ,&flag);
4920:     if (flag) {
4921:       MatMPIBAIJDiagonalScaleLocal(mat,diag);
4922:     } else {
4923:       int size;
4924:       MPI_Comm_size(mat->comm,&size);
4925:       if (size == 1) {
4926:         int n,m;
4927:         VecGetSize(diag,&n);
4928:         MatGetSize(mat,0,&m);
4929:         if (m == n) {
4930:           MatDiagonalScale(mat,0,diag);
4931:         } else {
4932:           SETERRQ(1,"Only supprted for sequential matrices when no ghost points/periodic conditions");
4933:         }
4934:       } else {
4935:         SETERRQ(1,"Only supported for MPIAIJ and MPIBAIJ parallel matrices");
4936:       }
4937:     }
4938:   }
4939:   PetscLogEventEnd(MAT_Scale,mat,0,0,0);
4940:   return(0);
4941: }

4943: #undef __FUNCT__  
4945: /*@ 
4946:    MatGetInertia - Gets the inertia from a factored matrix

4948:    Collective on Mat

4950:    Input Parameter:
4951: .  mat - the matrix

4953:    Output Parameters:
4954: +   nneg - number of negative eigenvalues
4955: .   nzero - number of zero eigenvalues
4956: -   npos - number of positive eigenvalues

4958:    Level: advanced

4960:    Notes: Matrix must have been factored by MatCholeskyFactor()


4963: @*/
4964: int MatGetInertia(Mat mat,int *nneg,int *nzero,int *npos)
4965: {
4966:   int        ierr;

4971:   if (!mat->factor)    SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
4972:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Numeric factor mat is not assembled");
4973:   if (!mat->ops->getinertia) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
4974:   (*mat->ops->getinertia)(mat,nneg,nzero,npos);
4975:   return(0);
4976: }

4978: /* ----------------------------------------------------------------*/
4979: #undef __FUNCT__  
4981: /*@
4982:    MatSolves - Solves A x = b, given a factored matrix, for a collection of vectors

4984:    Collective on Mat and Vecs

4986:    Input Parameters:
4987: +  mat - the factored matrix
4988: -  b - the right-hand-side vectors

4990:    Output Parameter:
4991: .  x - the result vectors

4993:    Notes:
4994:    The vectors b and x cannot be the same.  I.e., one cannot
4995:    call MatSolves(A,x,x).

4997:    Notes:
4998:    Most users should employ the simplified SLES interface for linear solvers
4999:    instead of working directly with matrix algebra routines such as this.
5000:    See, e.g., SLESCreate().

5002:    Level: developer

5004:    Concepts: matrices^triangular solves

5006: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd(), MatSolve()
5007: @*/
5008: int MatSolves(Mat mat,Vecs b,Vecs x)
5009: {

5015:   MatPreallocated(mat);
5016:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
5017:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
5018:   if (mat->M == 0 && mat->N == 0) return(0);

5020:   if (!mat->ops->solves) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
5021:   PetscLogEventBegin(MAT_Solves,mat,0,0,0);
5022:   (*mat->ops->solves)(mat,b,x);
5023:   PetscLogEventEnd(MAT_Solves,mat,0,0,0);
5024:   return(0);
5025: }