Actual source code: matrix.c

      1: #define PETSCMAT_DLL

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

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

     10: /* Logging support */
     11: PetscCookie PETSCMAT_DLLEXPORT MAT_COOKIE = 0;
     12: PetscEvent  MAT_Mult = 0, MAT_Mults = 0, MAT_MultConstrained = 0, MAT_MultAdd = 0, MAT_MultTranspose = 0;
     13: PetscEvent  MAT_MultTransposeConstrained = 0, MAT_MultTransposeAdd = 0, MAT_Solve = 0, MAT_Solves = 0, MAT_SolveAdd = 0, MAT_SolveTranspose = 0;
     14: PetscEvent  MAT_SolveTransposeAdd = 0, MAT_Relax = 0, MAT_ForwardSolve = 0, MAT_BackwardSolve = 0, MAT_LUFactor = 0, MAT_LUFactorSymbolic = 0;
     15: PetscEvent  MAT_LUFactorNumeric = 0, MAT_CholeskyFactor = 0, MAT_CholeskyFactorSymbolic = 0, MAT_CholeskyFactorNumeric = 0, MAT_ILUFactor = 0;
     16: PetscEvent  MAT_ILUFactorSymbolic = 0, MAT_ICCFactorSymbolic = 0, MAT_Copy = 0, MAT_Convert = 0, MAT_Scale = 0, MAT_AssemblyBegin = 0;
     17: PetscEvent  MAT_AssemblyEnd = 0, MAT_SetValues = 0, MAT_GetValues = 0, MAT_GetRow = 0, MAT_GetSubMatrices = 0, MAT_GetColoring = 0, MAT_GetOrdering = 0;
     18: PetscEvent  MAT_IncreaseOverlap = 0, MAT_Partitioning = 0, MAT_ZeroEntries = 0, MAT_Load = 0, MAT_View = 0, MAT_AXPY = 0, MAT_FDColoringCreate = 0;
     19: PetscEvent  MAT_FDColoringApply = 0,MAT_Transpose = 0,MAT_FDColoringFunction = 0;
     20: PetscEvent  MAT_MatMult = 0, MAT_MatMultSymbolic = 0, MAT_MatMultNumeric = 0;
     21: PetscEvent  MAT_PtAP = 0, MAT_PtAPSymbolic = 0, MAT_PtAPNumeric = 0;
     22: PetscEvent  MAT_MatMultTranspose = 0, MAT_MatMultTransposeSymbolic = 0, MAT_MatMultTransposeNumeric = 0;

     24: /* nasty global values for MatSetValue() */
     25: PetscInt    PETSCMAT_DLLEXPORT MatSetValue_Row = 0;
     26: PetscInt    PETSCMAT_DLLEXPORT MatSetValue_Column = 0;
     27: PetscScalar PETSCMAT_DLLEXPORT MatSetValue_Value = 0.0;

     31: /*@C
     32:    MatGetRow - Gets a row of a matrix.  You MUST call MatRestoreRow()
     33:    for each row that you get to ensure that your application does
     34:    not bleed memory.

     36:    Not Collective

     38:    Input Parameters:
     39: +  mat - the matrix
     40: -  row - the row to get

     42:    Output Parameters:
     43: +  ncols -  if not NULL, the number of nonzeros in the row
     44: .  cols - if not NULL, the column numbers
     45: -  vals - if not NULL, the values

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

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

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

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

     62:    You can only have one call to MatGetRow() outstanding for a particular
     63:    matrix at a time, per processor. MatGetRow() can only obtain rows
     64:    associated with the given processor, it cannot get rows from the
     65:    other processors; for that we suggest using MatGetSubMatrices(), then
     66:    MatGetRow() on the submatrix. The row indix passed to MatGetRows()
     67:    is in the global number of rows.

     69:    Fortran Notes:
     70:    The calling sequence from Fortran is
     71: .vb
     72:    MatGetRow(matrix,row,ncols,cols,values,ierr)
     73:          Mat     matrix (input)
     74:          integer row    (input)
     75:          integer ncols  (output)
     76:          integer cols(maxcols) (output)
     77:          double precision (or double complex) values(maxcols) output
     78: .ve
     79:    where maxcols >= maximum nonzeros in any row of the matrix.

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

     86:    Level: advanced

     88:    Concepts: matrices^row access

     90: .seealso: MatRestoreRow(), MatSetValues(), MatGetValues(), MatGetSubmatrices(), MatGetDiagonal()
     91: @*/

     93: PetscErrorCode PETSCMAT_DLLEXPORT MatGetRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
     94: {
     96:   PetscInt       incols;

    101:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    102:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    103:   if (!mat->ops->getrow) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    104:   MatPreallocated(mat);
    105:   PetscLogEventBegin(MAT_GetRow,mat,0,0,0);
    106:   (*mat->ops->getrow)(mat,row,&incols,(PetscInt **)cols,(PetscScalar **)vals);
    107:   if (ncols) *ncols = incols;
    108:   PetscLogEventEnd(MAT_GetRow,mat,0,0,0);
    109:   return(0);
    110: }

    114: /*@
    115:    MatConjugate - replaces the matrix values with their complex conjugates

    117:    Collective on Mat

    119:    Input Parameters:
    120: .  mat - the matrix

    122:    Level: advanced

    124: .seealso:  VecConjugate()
    125: @*/
    126: PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate(Mat mat)
    127: {

    132:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    133:   if (!mat->ops->conjugate) SETERRQ(PETSC_ERR_SUP,"Not provided for this matrix format, send email to petsc-maint@mcs.anl.gov");
    134:   (*mat->ops->conjugate)(mat);
    135:   return(0);
    136: }

    140: /*@C
    141:    MatRestoreRow - Frees any temporary space allocated by MatGetRow().

    143:    Not Collective

    145:    Input Parameters:
    146: +  mat - the matrix
    147: .  row - the row to get
    148: .  ncols, cols - the number of nonzeros and their columns
    149: -  vals - if nonzero the column values

    151:    Notes:
    152:    This routine should be called after you have finished examining the entries.

    154:    Fortran Notes:
    155:    The calling sequence from Fortran is
    156: .vb
    157:    MatRestoreRow(matrix,row,ncols,cols,values,ierr)
    158:       Mat     matrix (input)
    159:       integer row    (input)
    160:       integer ncols  (output)
    161:       integer cols(maxcols) (output)
    162:       double precision (or double complex) values(maxcols) output
    163: .ve
    164:    Where maxcols >= maximum nonzeros in any row of the matrix.

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

    169:    Level: advanced

    171: .seealso:  MatGetRow()
    172: @*/
    173: PetscErrorCode PETSCMAT_DLLEXPORT MatRestoreRow(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt *cols[],const PetscScalar *vals[])
    174: {

    180:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    181:   if (!mat->ops->restorerow) return(0);
    182:   (*mat->ops->restorerow)(mat,row,ncols,(PetscInt **)cols,(PetscScalar **)vals);
    183:   return(0);
    184: }

    188: /*@C
    189:    MatSetOptionsPrefix - Sets the prefix used for searching for all
    190:    Mat options in the database.

    192:    Collective on Mat

    194:    Input Parameter:
    195: +  A - the Mat context
    196: -  prefix - the prefix to prepend to all option names

    198:    Notes:
    199:    A hyphen (-) must NOT be given at the beginning of the prefix name.
    200:    The first character of all runtime options is AUTOMATICALLY the hyphen.

    202:    Level: advanced

    204: .keywords: Mat, set, options, prefix, database

    206: .seealso: MatSetFromOptions()
    207: @*/
    208: PetscErrorCode PETSCMAT_DLLEXPORT MatSetOptionsPrefix(Mat A,const char prefix[])
    209: {

    214:   PetscObjectSetOptionsPrefix((PetscObject)A,prefix);
    215:   return(0);
    216: }

    220: /*@C
    221:    MatAppendOptionsPrefix - Appends to the prefix used for searching for all
    222:    Mat options in the database.

    224:    Collective on Mat

    226:    Input Parameters:
    227: +  A - the Mat context
    228: -  prefix - the prefix to prepend to all option names

    230:    Notes:
    231:    A hyphen (-) must NOT be given at the beginning of the prefix name.
    232:    The first character of all runtime options is AUTOMATICALLY the hyphen.

    234:    Level: advanced

    236: .keywords: Mat, append, options, prefix, database

    238: .seealso: MatGetOptionsPrefix()
    239: @*/
    240: PetscErrorCode PETSCMAT_DLLEXPORT MatAppendOptionsPrefix(Mat A,const char prefix[])
    241: {
    243:
    246:   PetscObjectAppendOptionsPrefix((PetscObject)A,prefix);
    247:   return(0);
    248: }

    252: /*@C
    253:    MatGetOptionsPrefix - Sets the prefix used for searching for all
    254:    Mat options in the database.

    256:    Not Collective

    258:    Input Parameter:
    259: .  A - the Mat context

    261:    Output Parameter:
    262: .  prefix - pointer to the prefix string used

    264:    Notes: On the fortran side, the user should pass in a string 'prefix' of
    265:    sufficient length to hold the prefix.

    267:    Level: advanced

    269: .keywords: Mat, get, options, prefix, database

    271: .seealso: MatAppendOptionsPrefix()
    272: @*/
    273: PetscErrorCode PETSCMAT_DLLEXPORT MatGetOptionsPrefix(Mat A,const char *prefix[])
    274: {

    279:   PetscObjectGetOptionsPrefix((PetscObject)A,prefix);
    280:   return(0);
    281: }

    285: /*@
    286:    MatSetUp - Sets up the internal matrix data structures for the later use.

    288:    Collective on Mat

    290:    Input Parameters:
    291: .  A - the Mat context

    293:    Notes:
    294:    For basic use of the Mat classes the user need not explicitly call
    295:    MatSetUp(), since these actions will happen automatically.

    297:    Level: advanced

    299: .keywords: Mat, setup

    301: .seealso: MatCreate(), MatDestroy()
    302: @*/
    303: PetscErrorCode PETSCMAT_DLLEXPORT MatSetUp(Mat A)
    304: {

    309:   MatSetUpPreallocation(A);
    310:   MatSetFromOptions(A);
    311:   return(0);
    312: }

    316: /*@C
    317:    MatView - Visualizes a matrix object.

    319:    Collective on Mat

    321:    Input Parameters:
    322: +  mat - the matrix
    323: -  viewer - visualization context

    325:   Notes:
    326:   The available visualization contexts include
    327: +    PETSC_VIEWER_STDOUT_SELF - standard output (default)
    328: .    PETSC_VIEWER_STDOUT_WORLD - synchronized standard
    329:         output where only the first processor opens
    330:         the file.  All other processors send their
    331:         data to the first processor to print.
    332: -     PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure

    334:    The user can open alternative visualization contexts with
    335: +    PetscViewerASCIIOpen() - Outputs matrix to a specified file
    336: .    PetscViewerBinaryOpen() - Outputs matrix in binary to a
    337:          specified file; corresponding input uses MatLoad()
    338: .    PetscViewerDrawOpen() - Outputs nonzero matrix structure to
    339:          an X window display
    340: -    PetscViewerSocketOpen() - Outputs matrix to Socket viewer.
    341:          Currently only the sequential dense and AIJ
    342:          matrix types support the Socket viewer.

    344:    The user can call PetscViewerSetFormat() to specify the output
    345:    format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
    346:    PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen).  Available formats include
    347: +    PETSC_VIEWER_ASCII_DEFAULT - default, prints matrix contents
    348: .    PETSC_VIEWER_ASCII_MATLAB - prints matrix contents in Matlab format
    349: .    PETSC_VIEWER_ASCII_DENSE - prints entire matrix including zeros
    350: .    PETSC_VIEWER_ASCII_COMMON - prints matrix contents, using a sparse
    351:          format common among all matrix types
    352: .    PETSC_VIEWER_ASCII_IMPL - prints matrix contents, using an implementation-specific
    353:          format (which is in many cases the same as the default)
    354: .    PETSC_VIEWER_ASCII_INFO - prints basic information about the matrix
    355:          size and structure (not the matrix entries)
    356: .    PETSC_VIEWER_ASCII_INFO_DETAIL - prints more detailed information about
    357:          the matrix structure

    359:    Options Database Keys:
    360: +  -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
    361: .  -mat_view_info_detailed - Prints more detailed info
    362: .  -mat_view - Prints matrix in ASCII format
    363: .  -mat_view_matlab - Prints matrix in Matlab format
    364: .  -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
    365: .  -display <name> - Sets display name (default is host)
    366: .  -draw_pause <sec> - Sets number of seconds to pause after display
    367: .  -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
    368: .  -viewer_socket_machine <machine>
    369: .  -viewer_socket_port <port>
    370: .  -mat_view_binary - save matrix to file in binary format
    371: -  -viewer_binary_filename <name>
    372:    Level: beginner

    374:    Concepts: matrices^viewing
    375:    Concepts: matrices^plotting
    376:    Concepts: matrices^printing

    378: .seealso: PetscViewerSetFormat(), PetscViewerASCIIOpen(), PetscViewerDrawOpen(),
    379:           PetscViewerSocketOpen(), PetscViewerBinaryOpen(), MatLoad()
    380: @*/
    381: PetscErrorCode PETSCMAT_DLLEXPORT MatView(Mat mat,PetscViewer viewer)
    382: {
    383:   PetscErrorCode    ierr;
    384:   PetscInt          rows,cols;
    385:   PetscTruth        iascii;
    386:   const char        *cstr;
    387:   PetscViewerFormat format;

    392:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(mat->comm);
    395:   if (!mat->assembled) SETERRQ(PETSC_ERR_ORDER,"Must call MatAssemblyBegin/End() before viewing matrix");
    396:   MatPreallocated(mat);

    398:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    399:   if (iascii) {
    400:     PetscViewerGetFormat(viewer,&format);
    401:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
    402:       if (mat->prefix) {
    403:         PetscViewerASCIIPrintf(viewer,"Matrix Object:(%s)\n",mat->prefix);
    404:       } else {
    405:         PetscViewerASCIIPrintf(viewer,"Matrix Object:\n");
    406:       }
    407:       PetscViewerASCIIPushTab(viewer);
    408:       MatGetType(mat,&cstr);
    409:       MatGetSize(mat,&rows,&cols);
    410:       PetscViewerASCIIPrintf(viewer,"type=%s, rows=%D, cols=%D\n",cstr,rows,cols);
    411:       if (mat->ops->getinfo) {
    412:         MatInfo info;
    413:         MatGetInfo(mat,MAT_GLOBAL_SUM,&info);
    414:         PetscViewerASCIIPrintf(viewer,"total: nonzeros=%D, allocated nonzeros=%D\n",
    415:                           (PetscInt)info.nz_used,(PetscInt)info.nz_allocated);
    416:       }
    417:     }
    418:   }
    419:   if (mat->ops->view) {
    420:     PetscViewerASCIIPushTab(viewer);
    421:     (*mat->ops->view)(mat,viewer);
    422:     PetscViewerASCIIPopTab(viewer);
    423:   } else if (!iascii) {
    424:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
    425:   }
    426:   if (iascii) {
    427:     PetscViewerGetFormat(viewer,&format);
    428:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
    429:       PetscViewerASCIIPopTab(viewer);
    430:     }
    431:   }
    432:   return(0);
    433: }

    437: /*@C
    438:    MatScaleSystem - Scale a vector solution and right hand side to
    439:    match the scaling of a scaled matrix.
    440:
    441:    Collective on Mat

    443:    Input Parameter:
    444: +  mat - the matrix
    445: .  x - solution vector (or PETSC_NULL)
    446: -  b - right hand side vector (or PETSC_NULL)

    449:    Notes:
    450:    For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
    451:    internally scaled, so this does nothing. For MPIROWBS it
    452:    permutes and diagonally scales.

    454:    The KSP methods automatically call this routine when required
    455:    (via PCPreSolve()) so it is rarely used directly.

    457:    Level: Developer

    459:    Concepts: matrices^scaling

    461: .seealso: MatUseScaledForm(), MatUnScaleSystem()
    462: @*/
    463: PetscErrorCode PETSCMAT_DLLEXPORT MatScaleSystem(Mat mat,Vec x,Vec b)
    464: {

    470:   MatPreallocated(mat);

    474:   if (mat->ops->scalesystem) {
    475:     (*mat->ops->scalesystem)(mat,x,b);
    476:   }
    477:   return(0);
    478: }

    482: /*@C
    483:    MatUnScaleSystem - Unscales a vector solution and right hand side to
    484:    match the original scaling of a scaled matrix.
    485:
    486:    Collective on Mat

    488:    Input Parameter:
    489: +  mat - the matrix
    490: .  x - solution vector (or PETSC_NULL)
    491: -  b - right hand side vector (or PETSC_NULL)

    494:    Notes:
    495:    For AIJ, BAIJ, and BDiag matrix formats, the matrices are not
    496:    internally scaled, so this does nothing. For MPIROWBS it
    497:    permutes and diagonally scales.

    499:    The KSP methods automatically call this routine when required
    500:    (via PCPreSolve()) so it is rarely used directly.

    502:    Level: Developer

    504: .seealso: MatUseScaledForm(), MatScaleSystem()
    505: @*/
    506: PetscErrorCode PETSCMAT_DLLEXPORT MatUnScaleSystem(Mat mat,Vec x,Vec b)
    507: {

    513:   MatPreallocated(mat);
    516:   if (mat->ops->unscalesystem) {
    517:     (*mat->ops->unscalesystem)(mat,x,b);
    518:   }
    519:   return(0);
    520: }

    524: /*@C
    525:    MatUseScaledForm - For matrix storage formats that scale the
    526:    matrix (for example MPIRowBS matrices are diagonally scaled on
    527:    assembly) indicates matrix operations (MatMult() etc) are
    528:    applied using the scaled matrix.
    529:
    530:    Collective on Mat

    532:    Input Parameter:
    533: +  mat - the matrix
    534: -  scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for
    535:             applying the original matrix

    537:    Notes:
    538:    For scaled matrix formats, applying the original, unscaled matrix
    539:    will be slightly more expensive

    541:    Level: Developer

    543: .seealso: MatScaleSystem(), MatUnScaleSystem()
    544: @*/
    545: PetscErrorCode PETSCMAT_DLLEXPORT MatUseScaledForm(Mat mat,PetscTruth scaled)
    546: {

    552:   MatPreallocated(mat);
    553:   if (mat->ops->usescaledform) {
    554:     (*mat->ops->usescaledform)(mat,scaled);
    555:   }
    556:   return(0);
    557: }

    561: /*@C
    562:    MatDestroy - Frees space taken by a matrix.
    563:
    564:    Collective on Mat

    566:    Input Parameter:
    567: .  A - the matrix

    569:    Level: beginner

    571: @*/
    572: PetscErrorCode PETSCMAT_DLLEXPORT MatDestroy(Mat A)
    573: {

    578:   if (--A->refct > 0) return(0);

    581:   MatPreallocated(A);
    582:   /* if memory was published with AMS then destroy it */
    583:   PetscObjectDepublish(A);
    584:   if (A->mapping) {
    585:     ISLocalToGlobalMappingDestroy(A->mapping);
    586:   }
    587:   if (A->bmapping) {
    588:     ISLocalToGlobalMappingDestroy(A->bmapping);
    589:   }
    590:   if (A->rmap) {
    591:     PetscMapDestroy(A->rmap);
    592:   }
    593:   if (A->cmap) {
    594:     PetscMapDestroy(A->cmap);
    595:   }
    596:   (*A->ops->destroy)(A);
    597:   PetscHeaderDestroy(A);
    598:   return(0);
    599: }

    603: /*@
    604:    MatValid - Checks whether a matrix object is valid.

    606:    Collective on Mat

    608:    Input Parameter:
    609: .  m - the matrix to check

    611:    Output Parameter:
    612:    flg - flag indicating matrix status, either
    613:    PETSC_TRUE if matrix is valid, or PETSC_FALSE otherwise.

    615:    Level: developer

    617:    Concepts: matrices^validity
    618: @*/
    619: PetscErrorCode PETSCMAT_DLLEXPORT MatValid(Mat m,PetscTruth *flg)
    620: {
    623:   if (!m)                           *flg = PETSC_FALSE;
    624:   else if (m->cookie != MAT_COOKIE) *flg = PETSC_FALSE;
    625:   else                              *flg = PETSC_TRUE;
    626:   return(0);
    627: }

    631: /*@
    632:    MatSetValues - Inserts or adds a block of values into a matrix.
    633:    These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
    634:    MUST be called after all calls to MatSetValues() have been completed.

    636:    Not Collective

    638:    Input Parameters:
    639: +  mat - the matrix
    640: .  v - a logically two-dimensional array of values
    641: .  m, idxm - the number of rows and their global indices
    642: .  n, idxn - the number of columns and their global indices
    643: -  addv - either ADD_VALUES or INSERT_VALUES, where
    644:    ADD_VALUES adds values to any existing entries, and
    645:    INSERT_VALUES replaces existing entries with new values

    647:    Notes:
    648:    By default the values, v, are row-oriented and unsorted.
    649:    See MatSetOption() for other options.

    651:    Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES
    652:    options cannot be mixed without intervening calls to the assembly
    653:    routines.

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

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

    663:    Efficiency Alert:
    664:    The routine MatSetValuesBlocked() may offer much better efficiency
    665:    for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).

    667:    Level: beginner

    669:    Concepts: matrices^putting entries in

    671: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
    672:           InsertMode, INSERT_VALUES, ADD_VALUES
    673: @*/
    674: PetscErrorCode PETSCMAT_DLLEXPORT MatSetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
    675: {

    679:   if (!m || !n) return(0); /* no values to insert */
    685:   MatPreallocated(mat);
    686:   if (mat->insertmode == NOT_SET_VALUES) {
    687:     mat->insertmode = addv;
    688:   }
    689: #if defined(PETSC_USE_DEBUG)
    690:   else if (mat->insertmode != addv) {
    691:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
    692:   }
    693:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    694: #endif

    696:   if (mat->assembled) {
    697:     mat->was_assembled = PETSC_TRUE;
    698:     mat->assembled     = PETSC_FALSE;
    699:   }
    700:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
    701:   if (!mat->ops->setvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    702:   (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);
    703:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
    704:   return(0);
    705: }

    709: /*@C
    710:    MatSetValuesStencil - Inserts or adds a block of values into a matrix.
    711:      Using structured grid indexing

    713:    Not Collective

    715:    Input Parameters:
    716: +  mat - the matrix
    717: .  v - a logically two-dimensional array of values
    718: .  m - number of rows being entered
    719: .  idxm - grid coordinates (and component number when dof > 1) for matrix rows being entered
    720: .  n - number of columns being entered
    721: .  idxn - grid coordinates (and component number when dof > 1) for matrix columns being entered
    722: -  addv - either ADD_VALUES or INSERT_VALUES, where
    723:    ADD_VALUES adds values to any existing entries, and
    724:    INSERT_VALUES replaces existing entries with new values

    726:    Notes:
    727:    By default the values, v, are row-oriented and unsorted.
    728:    See MatSetOption() for other options.

    730:    Calls to MatSetValuesStencil() with the INSERT_VALUES and ADD_VALUES
    731:    options cannot be mixed without intervening calls to the assembly
    732:    routines.

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

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

    739:    For setting/accessing vector values via array coordinates you can use the DAVecGetArray() routine

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

    744:    The columns and rows in the stencil passed in MUST be contained within the
    745:    ghost region of the given process as set with DACreateXXX() or MatSetStencil(). For example,
    746:    if you create a DA with an overlap of one grid level and on a particular process its first
    747:    local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
    748:    first i index you can use in your column and row indices in MatSetStencil() is 5.

    750:    In Fortran idxm and idxn should be declared as
    751: $     MatStencil idxm(4,m),idxn(4,n)
    752:    and the values inserted using
    753: $    idxm(MatStencil_i,1) = i
    754: $    idxm(MatStencil_j,1) = j
    755: $    idxm(MatStencil_k,1) = k
    756: $    idxm(MatStencil_c,1) = c
    757:    etc

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

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

    767:    Efficiency Alert:
    768:    The routine MatSetValuesBlockedStencil() may offer much better efficiency
    769:    for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).

    771:    Level: beginner

    773:    Concepts: matrices^putting entries in

    775: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
    776:           MatSetValues(), MatSetValuesBlockedStencil(), MatSetStencil(), DAGetMatrix(), DAVecGetArray(), MatStencil
    777: @*/
    778: PetscErrorCode PETSCMAT_DLLEXPORT MatSetValuesStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
    779: {
    781:   PetscInt       j,i,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
    782:   PetscInt       *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);

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

    792:   if (m > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 128 rows at a time; trying to set %D",m);
    793:   if (n > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 256 columns at a time; trying to set %D",n);

    795:   for (i=0; i<m; i++) {
    796:     for (j=0; j<3-sdim; j++) dxm++;
    797:     tmp = *dxm++ - starts[0];
    798:     for (j=0; j<dim-1; j++) {
    799:       if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
    800:       else                                       tmp = tmp*dims[j] + dxm[-1] - starts[j+1];
    801:     }
    802:     if (mat->stencil.noc) dxm++;
    803:     jdxm[i] = tmp;
    804:   }
    805:   for (i=0; i<n; i++) {
    806:     for (j=0; j<3-sdim; j++) dxn++;
    807:     tmp = *dxn++ - starts[0];
    808:     for (j=0; j<dim-1; j++) {
    809:       if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
    810:       else                                       tmp = tmp*dims[j] + dxn[-1] - starts[j+1];
    811:     }
    812:     if (mat->stencil.noc) dxn++;
    813:     jdxn[i] = tmp;
    814:   }
    815:   MatSetValuesLocal(mat,m,jdxm,n,jdxn,v,addv);
    816:   return(0);
    817: }

    821: /*@C
    822:    MatSetValuesBlockedStencil - Inserts or adds a block of values into a matrix.
    823:      Using structured grid indexing

    825:    Not Collective

    827:    Input Parameters:
    828: +  mat - the matrix
    829: .  v - a logically two-dimensional array of values
    830: .  m - number of rows being entered
    831: .  idxm - grid coordinates for matrix rows being entered
    832: .  n - number of columns being entered
    833: .  idxn - grid coordinates for matrix columns being entered
    834: -  addv - either ADD_VALUES or INSERT_VALUES, where
    835:    ADD_VALUES adds values to any existing entries, and
    836:    INSERT_VALUES replaces existing entries with new values

    838:    Notes:
    839:    By default the values, v, are row-oriented and unsorted.
    840:    See MatSetOption() for other options.

    842:    Calls to MatSetValuesBlockedStencil() with the INSERT_VALUES and ADD_VALUES
    843:    options cannot be mixed without intervening calls to the assembly
    844:    routines.

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

    848:    MatSetValuesBlockedStencil() uses 0-based row and column numbers in Fortran
    849:    as well as in C.

    851:    For setting/accessing vector values via array coordinates you can use the DAVecGetArray() routine

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

    856:    The columns and rows in the stencil passed in MUST be contained within the
    857:    ghost region of the given process as set with DACreateXXX() or MatSetStencil(). For example,
    858:    if you create a DA with an overlap of one grid level and on a particular process its first
    859:    local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the
    860:    first i index you can use in your column and row indices in MatSetStencil() is 5.

    862:    In Fortran idxm and idxn should be declared as
    863: $     MatStencil idxm(4,m),idxn(4,n)
    864:    and the values inserted using
    865: $    idxm(MatStencil_i,1) = i
    866: $    idxm(MatStencil_j,1) = j
    867: $    idxm(MatStencil_k,1) = k
    868:    etc

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

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

    878:    Level: beginner

    880:    Concepts: matrices^putting entries in

    882: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
    883:           MatSetValues(), MatSetValuesStencil(), MatSetStencil(), DAGetMatrix(), DAVecGetArray(), MatStencil
    884: @*/
    885: PetscErrorCode PETSCMAT_DLLEXPORT MatSetValuesBlockedStencil(Mat mat,PetscInt m,const MatStencil idxm[],PetscInt n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv)
    886: {
    888:   PetscInt       j,i,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp;
    889:   PetscInt       *starts = mat->stencil.starts,*dxm = (PetscInt*)idxm,*dxn = (PetscInt*)idxn,sdim = dim - (1 - (PetscInt)mat->stencil.noc);

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

    899:   if (m > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 128 rows at a time; trying to set %D",m);
    900:   if (n > 128) SETERRQ1(PETSC_ERR_SUP,"Can only set 256 columns at a time; trying to set %D",n);

    902:   for (i=0; i<m; i++) {
    903:     for (j=0; j<3-sdim; j++) dxm++;
    904:     tmp = *dxm++ - starts[0];
    905:     for (j=0; j<sdim-1; j++) {
    906:       if ((*dxm++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
    907:       else                                      tmp = tmp*dims[j] + dxm[-1] - starts[j+1];
    908:     }
    909:     dxm++;
    910:     jdxm[i] = tmp;
    911:   }
    912:   for (i=0; i<n; i++) {
    913:     for (j=0; j<3-sdim; j++) dxn++;
    914:     tmp = *dxn++ - starts[0];
    915:     for (j=0; j<sdim-1; j++) {
    916:       if ((*dxn++ - starts[j+1]) < 0 || tmp < 0) tmp = PETSC_MIN_INT;
    917:       else                                       tmp = tmp*dims[j] + dxn[-1] - starts[j+1];
    918:     }
    919:     dxn++;
    920:     jdxn[i] = tmp;
    921:   }
    922:   MatSetValuesBlockedLocal(mat,m,jdxm,n,jdxn,v,addv);
    923:   return(0);
    924: }

    928: /*@
    929:    MatSetStencil - Sets the grid information for setting values into a matrix via
    930:         MatSetValuesStencil()

    932:    Not Collective

    934:    Input Parameters:
    935: +  mat - the matrix
    936: .  dim - dimension of the grid 1, 2, or 3
    937: .  dims - number of grid points in x, y, and z direction, including ghost points on your processor
    938: .  starts - starting point of ghost nodes on your processor in x, y, and z direction
    939: -  dof - number of degrees of freedom per node

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

    945:    For matrices generated with DAGetMatrix() this routine is automatically called and so not needed by the
    946:    user.
    947:
    948:    Level: beginner

    950:    Concepts: matrices^putting entries in

    952: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal()
    953:           MatSetValues(), MatSetValuesBlockedStencil(), MatSetValuesStencil()
    954: @*/
    955: PetscErrorCode PETSCMAT_DLLEXPORT MatSetStencil(Mat mat,PetscInt dim,const PetscInt dims[],const PetscInt starts[],PetscInt dof)
    956: {
    957:   PetscInt i;

    964:   mat->stencil.dim = dim + (dof > 1);
    965:   for (i=0; i<dim; i++) {
    966:     mat->stencil.dims[i]   = dims[dim-i-1];      /* copy the values in backwards */
    967:     mat->stencil.starts[i] = starts[dim-i-1];
    968:   }
    969:   mat->stencil.dims[dim]   = dof;
    970:   mat->stencil.starts[dim] = 0;
    971:   mat->stencil.noc         = (PetscTruth)(dof == 1);
    972:   return(0);
    973: }

    977: /*@
    978:    MatSetValuesBlocked - Inserts or adds a block of values into a matrix.

    980:    Not Collective

    982:    Input Parameters:
    983: +  mat - the matrix
    984: .  v - a logically two-dimensional array of values
    985: .  m, idxm - the number of block rows and their global block indices
    986: .  n, idxn - the number of block columns and their global block indices
    987: -  addv - either ADD_VALUES or INSERT_VALUES, where
    988:    ADD_VALUES adds values to any existing entries, and
    989:    INSERT_VALUES replaces existing entries with new values

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

    995:    Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
    996:    options cannot be mixed without intervening calls to the assembly
    997:    routines.

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

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

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

    1013:    Restrictions:
    1014:    MatSetValuesBlocked() is currently supported only for the BAIJ and SBAIJ formats

    1016:    Level: intermediate

    1018:    Concepts: matrices^putting entries in blocked

    1020: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal()
    1021: @*/
    1022: PetscErrorCode PETSCMAT_DLLEXPORT MatSetValuesBlocked(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
    1023: {

    1027:   if (!m || !n) return(0); /* no values to insert */
    1033:   MatPreallocated(mat);
    1034:   if (mat->insertmode == NOT_SET_VALUES) {
    1035:     mat->insertmode = addv;
    1036:   }
    1037: #if defined(PETSC_USE_DEBUG)
    1038:   else if (mat->insertmode != addv) {
    1039:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
    1040:   }
    1041:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1042: #endif

    1044:   if (mat->assembled) {
    1045:     mat->was_assembled = PETSC_TRUE;
    1046:     mat->assembled     = PETSC_FALSE;
    1047:   }
    1048:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
    1049:   if (!mat->ops->setvaluesblocked) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    1050:   (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);
    1051:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
    1052:   return(0);
    1053: }

    1057: /*@
    1058:    MatGetValues - Gets a block of values from a matrix.

    1060:    Not Collective; currently only returns a local block

    1062:    Input Parameters:
    1063: +  mat - the matrix
    1064: .  v - a logically two-dimensional array for storing the values
    1065: .  m, idxm - the number of rows and their global indices
    1066: -  n, idxn - the number of columns and their global indices

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

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

    1076:    MatGetValues() requires that the matrix has been assembled
    1077:    with MatAssemblyBegin()/MatAssemblyEnd().  Thus, calls to
    1078:    MatSetValues() and MatGetValues() CANNOT be made in succession
    1079:    without intermediate matrix assembly.

    1081:    Level: advanced

    1083:    Concepts: matrices^accessing values

    1085: .seealso: MatGetRow(), MatGetSubmatrices(), MatSetValues()
    1086: @*/
    1087: PetscErrorCode PETSCMAT_DLLEXPORT MatGetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
    1088: {

    1097:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    1098:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1099:   if (!mat->ops->getvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    1100:   MatPreallocated(mat);

    1102:   PetscLogEventBegin(MAT_GetValues,mat,0,0,0);
    1103:   (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);
    1104:   PetscLogEventEnd(MAT_GetValues,mat,0,0,0);
    1105:   return(0);
    1106: }

    1110: /*@
    1111:    MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by
    1112:    the routine MatSetValuesLocal() to allow users to insert matrix entries
    1113:    using a local (per-processor) numbering.

    1115:    Not Collective

    1117:    Input Parameters:
    1118: +  x - the matrix
    1119: -  mapping - mapping created with ISLocalToGlobalMappingCreate()
    1120:              or ISLocalToGlobalMappingCreateIS()

    1122:    Level: intermediate

    1124:    Concepts: matrices^local to global mapping
    1125:    Concepts: local to global mapping^for matrices

    1127: .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal()
    1128: @*/
    1129: PetscErrorCode PETSCMAT_DLLEXPORT MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping)
    1130: {
    1136:   if (x->mapping) {
    1137:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
    1138:   }
    1139:   MatPreallocated(x);

    1141:   if (x->ops->setlocaltoglobalmapping) {
    1142:     (*x->ops->setlocaltoglobalmapping)(x,mapping);
    1143:   } else {
    1144:     x->mapping = mapping;
    1145:     PetscObjectReference((PetscObject)mapping);
    1146:   }
    1147:   return(0);
    1148: }

    1152: /*@
    1153:    MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use
    1154:    by the routine MatSetValuesBlockedLocal() to allow users to insert matrix
    1155:    entries using a local (per-processor) numbering.

    1157:    Not Collective

    1159:    Input Parameters:
    1160: +  x - the matrix
    1161: -  mapping - mapping created with ISLocalToGlobalMappingCreate() or
    1162:              ISLocalToGlobalMappingCreateIS()

    1164:    Level: intermediate

    1166:    Concepts: matrices^local to global mapping blocked
    1167:    Concepts: local to global mapping^for matrices, blocked

    1169: .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(),
    1170:            MatSetValuesBlocked(), MatSetValuesLocal()
    1171: @*/
    1172: PetscErrorCode PETSCMAT_DLLEXPORT MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping mapping)
    1173: {
    1179:   if (x->bmapping) {
    1180:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix");
    1181:   }
    1182:   x->bmapping = mapping;
    1183:   PetscObjectReference((PetscObject)mapping);
    1184:   return(0);
    1185: }

    1189: /*@
    1190:    MatSetValuesLocal - Inserts or adds values into certain locations of a matrix,
    1191:    using a local ordering of the nodes.

    1193:    Not Collective

    1195:    Input Parameters:
    1196: +  x - the matrix
    1197: .  nrow, irow - number of rows and their local indices
    1198: .  ncol, icol - number of columns and their local indices
    1199: .  y -  a logically two-dimensional array of values
    1200: -  addv - either INSERT_VALUES or ADD_VALUES, where
    1201:    ADD_VALUES adds values to any existing entries, and
    1202:    INSERT_VALUES replaces existing entries with new values

    1204:    Notes:
    1205:    Before calling MatSetValuesLocal(), the user must first set the
    1206:    local-to-global mapping by calling MatSetLocalToGlobalMapping().

    1208:    Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES
    1209:    options cannot be mixed without intervening calls to the assembly
    1210:    routines.

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

    1215:    Level: intermediate

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

    1219: .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping(),
    1220:            MatSetValueLocal()
    1221: @*/
    1222: PetscErrorCode PETSCMAT_DLLEXPORT MatSetValuesLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
    1223: {
    1225:   PetscInt       irowm[2048],icolm[2048];

    1233:   MatPreallocated(mat);
    1234:   if (mat->insertmode == NOT_SET_VALUES) {
    1235:     mat->insertmode = addv;
    1236:   }
    1237: #if defined(PETSC_USE_DEBUG)
    1238:   else if (mat->insertmode != addv) {
    1239:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
    1240:   }
    1241:   if (!mat->ops->setvalueslocal && (nrow > 2048 || ncol > 2048)) {
    1242:     SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %D %D",nrow,ncol);
    1243:   }
    1244:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1245: #endif

    1247:   if (mat->assembled) {
    1248:     mat->was_assembled = PETSC_TRUE;
    1249:     mat->assembled     = PETSC_FALSE;
    1250:   }
    1251:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
    1252:   if (!mat->ops->setvalueslocal) {
    1253:     ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm);
    1254:     ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);
    1255:     (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);
    1256:   } else {
    1257:     (*mat->ops->setvalueslocal)(mat,nrow,irow,ncol,icol,y,addv);
    1258:   }
    1259:   mat->same_nonzero = PETSC_FALSE;
    1260:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
    1261:   return(0);
    1262: }

    1266: /*@
    1267:    MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix,
    1268:    using a local ordering of the nodes a block at a time.

    1270:    Not Collective

    1272:    Input Parameters:
    1273: +  x - the matrix
    1274: .  nrow, irow - number of rows and their local indices
    1275: .  ncol, icol - number of columns and their local indices
    1276: .  y -  a logically two-dimensional array of values
    1277: -  addv - either INSERT_VALUES or ADD_VALUES, where
    1278:    ADD_VALUES adds values to any existing entries, and
    1279:    INSERT_VALUES replaces existing entries with new values

    1281:    Notes:
    1282:    Before calling MatSetValuesBlockedLocal(), the user must first set the
    1283:    local-to-global mapping by calling MatSetLocalToGlobalMappingBlock(),
    1284:    where the mapping MUST be set for matrix blocks, not for matrix elements.

    1286:    Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
    1287:    options cannot be mixed without intervening calls to the assembly
    1288:    routines.

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

    1293:    Level: intermediate

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

    1297: .seealso:  MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked()
    1298: @*/
    1299: PetscErrorCode PETSCMAT_DLLEXPORT MatSetValuesBlockedLocal(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
    1300: {
    1302:   PetscInt       irowm[2048],icolm[2048];

    1310:   MatPreallocated(mat);
    1311:   if (mat->insertmode == NOT_SET_VALUES) {
    1312:     mat->insertmode = addv;
    1313:   }
    1314: #if defined(PETSC_USE_DEBUG)
    1315:   else if (mat->insertmode != addv) {
    1316:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values");
    1317:   }
    1318:   if (!mat->bmapping) {
    1319:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with MatSetLocalToGlobalMappingBlock()");
    1320:   }
    1321:   if (nrow > 2048 || ncol > 2048) {
    1322:     SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %D %D",nrow,ncol);
    1323:   }
    1324:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1325: #endif

    1327:   if (mat->assembled) {
    1328:     mat->was_assembled = PETSC_TRUE;
    1329:     mat->assembled     = PETSC_FALSE;
    1330:   }
    1331:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
    1332:   ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm);
    1333:   ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm);
    1334:   (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);
    1335:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
    1336:   return(0);
    1337: }

    1339: /* --------------------------------------------------------*/
    1342: /*@
    1343:    MatMult - Computes the matrix-vector product, y = Ax.

    1345:    Collective on Mat and Vec

    1347:    Input Parameters:
    1348: +  mat - the matrix
    1349: -  x   - the vector to be multiplied

    1351:    Output Parameters:
    1352: .  y - the result

    1354:    Notes:
    1355:    The vectors x and y cannot be the same.  I.e., one cannot
    1356:    call MatMult(A,y,y).

    1358:    Level: beginner

    1360:    Concepts: matrix-vector product

    1362: .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd()
    1363: @*/
    1364: PetscErrorCode PETSCMAT_DLLEXPORT MatMult(Mat mat,Vec x,Vec y)
    1365: {

    1374:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    1375:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1376:   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
    1377: #ifndef PETSC_HAVE_CONSTRAINTS
    1378:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->N,x->N);
    1379:   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->M,y->N);
    1380:   if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %D %D",mat->m,y->n);
    1381: #endif
    1382:   MatPreallocated(mat);

    1384:   if (mat->nullsp) {
    1385:     MatNullSpaceRemove(mat->nullsp,x,&x);
    1386:   }

    1388:   PetscLogEventBegin(MAT_Mult,mat,x,y,0);
    1389:   (*mat->ops->mult)(mat,x,y);
    1390:   PetscLogEventEnd(MAT_Mult,mat,x,y,0);

    1392:   if (mat->nullsp) {
    1393:     MatNullSpaceRemove(mat->nullsp,y,PETSC_NULL);
    1394:   }
    1395:   PetscObjectStateIncrease((PetscObject)y);
    1396:   return(0);
    1397: }

    1401: /*@
    1402:    MatMultTranspose - Computes matrix transpose times a vector.

    1404:    Collective on Mat and Vec

    1406:    Input Parameters:
    1407: +  mat - the matrix
    1408: -  x   - the vector to be multilplied

    1410:    Output Parameters:
    1411: .  y - the result

    1413:    Notes:
    1414:    The vectors x and y cannot be the same.  I.e., one cannot
    1415:    call MatMultTranspose(A,y,y).

    1417:    Level: beginner

    1419:    Concepts: matrix vector product^transpose

    1421: .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd()
    1422: @*/
    1423: PetscErrorCode PETSCMAT_DLLEXPORT MatMultTranspose(Mat mat,Vec x,Vec y)
    1424: {
    1426:   PetscTruth     flg1, flg2;

    1434:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    1435:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1436:   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
    1437: #ifndef PETSC_HAVE_CONSTRAINTS
    1438:   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->M,x->N);
    1439:   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->N,y->N);
    1440: #endif
    1441:   MatPreallocated(mat);

    1443:   if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP, "Operation not supported");
    1444:   PetscLogEventBegin(MAT_MultTranspose,mat,x,y,0);
    1445:   if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply tranpose defined");
    1446:
    1447:   PetscTypeCompare((PetscObject)mat,MATSEQSBAIJ,&flg1);
    1448:   PetscTypeCompare((PetscObject)mat,MATMPISBAIJ,&flg2);
    1449:   if (flg1 || flg2) { /* mat is in sbaij format */
    1450:     (*mat->ops->mult)(mat,x,y);
    1451:   } else {
    1452:     (*mat->ops->multtranspose)(mat,x,y);
    1453:   }
    1454:   PetscLogEventEnd(MAT_MultTranspose,mat,x,y,0);
    1455:   PetscObjectStateIncrease((PetscObject)y);
    1456:   return(0);
    1457: }

    1461: /*@
    1462:     MatMultAdd -  Computes v3 = v2 + A * v1.

    1464:     Collective on Mat and Vec

    1466:     Input Parameters:
    1467: +   mat - the matrix
    1468: -   v1, v2 - the vectors

    1470:     Output Parameters:
    1471: .   v3 - the result

    1473:     Notes:
    1474:     The vectors v1 and v3 cannot be the same.  I.e., one cannot
    1475:     call MatMultAdd(A,v1,v2,v1).

    1477:     Level: beginner

    1479:     Concepts: matrix vector product^addition

    1481: .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd()
    1482: @*/
    1483: PetscErrorCode PETSCMAT_DLLEXPORT MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3)
    1484: {

    1494:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    1495:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1496:   if (mat->N != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %D %D",mat->N,v1->N);
    1497:   if (mat->M != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %D %D",mat->M,v2->N);
    1498:   if (mat->M != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %D %D",mat->M,v3->N);
    1499:   if (mat->m != v3->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: local dim %D %D",mat->m,v3->n);
    1500:   if (mat->m != v2->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: local dim %D %D",mat->m,v2->n);
    1501:   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
    1502:   MatPreallocated(mat);

    1504:   PetscLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);
    1505:   (*mat->ops->multadd)(mat,v1,v2,v3);
    1506:   PetscLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);
    1507:   PetscObjectStateIncrease((PetscObject)v3);
    1508:   return(0);
    1509: }

    1513: /*@
    1514:    MatMultTransposeAdd - Computes v3 = v2 + A' * v1.

    1516:    Collective on Mat and Vec

    1518:    Input Parameters:
    1519: +  mat - the matrix
    1520: -  v1, v2 - the vectors

    1522:    Output Parameters:
    1523: .  v3 - the result

    1525:    Notes:
    1526:    The vectors v1 and v3 cannot be the same.  I.e., one cannot
    1527:    call MatMultTransposeAdd(A,v1,v2,v1).

    1529:    Level: beginner

    1531:    Concepts: matrix vector product^transpose and addition

    1533: .seealso: MatMultTranspose(), MatMultAdd(), MatMult()
    1534: @*/
    1535: PetscErrorCode PETSCMAT_DLLEXPORT MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3)
    1536: {

    1546:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    1547:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1548:   if (!mat->ops->multtransposeadd) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    1549:   if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors");
    1550:   if (mat->M != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %D %D",mat->M,v1->N);
    1551:   if (mat->N != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %D %D",mat->N,v2->N);
    1552:   if (mat->N != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %D %D",mat->N,v3->N);
    1553:   MatPreallocated(mat);

    1555:   PetscLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);
    1556:   (*mat->ops->multtransposeadd)(mat,v1,v2,v3);
    1557:   PetscLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);
    1558:   PetscObjectStateIncrease((PetscObject)v3);
    1559:   return(0);
    1560: }

    1564: /*@
    1565:    MatMultConstrained - The inner multiplication routine for a
    1566:    constrained matrix P^T A P.

    1568:    Collective on Mat and Vec

    1570:    Input Parameters:
    1571: +  mat - the matrix
    1572: -  x   - the vector to be multilplied

    1574:    Output Parameters:
    1575: .  y - the result

    1577:    Notes:
    1578:    The vectors x and y cannot be the same.  I.e., one cannot
    1579:    call MatMult(A,y,y).

    1581:    Level: beginner

    1583: .keywords: matrix, multiply, matrix-vector product, constraint
    1584: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
    1585: @*/
    1586: PetscErrorCode PETSCMAT_DLLEXPORT MatMultConstrained(Mat mat,Vec x,Vec y)
    1587: {

    1594:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    1595:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1596:   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
    1597:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->N,x->N);
    1598:   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->M,y->N);
    1599:   if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %D %D",mat->m,y->n);

    1601:   PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
    1602:   (*mat->ops->multconstrained)(mat,x,y);
    1603:   PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
    1604:   PetscObjectStateIncrease((PetscObject)y);

    1606:   return(0);
    1607: }

    1611: /*@
    1612:    MatMultTransposeConstrained - The inner multiplication routine for a
    1613:    constrained matrix P^T A^T P.

    1615:    Collective on Mat and Vec

    1617:    Input Parameters:
    1618: +  mat - the matrix
    1619: -  x   - the vector to be multilplied

    1621:    Output Parameters:
    1622: .  y - the result

    1624:    Notes:
    1625:    The vectors x and y cannot be the same.  I.e., one cannot
    1626:    call MatMult(A,y,y).

    1628:    Level: beginner

    1630: .keywords: matrix, multiply, matrix-vector product, constraint
    1631: .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd()
    1632: @*/
    1633: PetscErrorCode PETSCMAT_DLLEXPORT MatMultTransposeConstrained(Mat mat,Vec x,Vec y)
    1634: {

    1641:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    1642:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1643:   if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors");
    1644:   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->N,x->N);
    1645:   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->M,y->N);

    1647:   PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);
    1648:   (*mat->ops->multtransposeconstrained)(mat,x,y);
    1649:   PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);
    1650:   PetscObjectStateIncrease((PetscObject)y);

    1652:   return(0);
    1653: }
    1654: /* ------------------------------------------------------------*/
    1657: /*@C
    1658:    MatGetInfo - Returns information about matrix storage (number of
    1659:    nonzeros, memory, etc.).

    1661:    Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used
    1662:    as the flag

    1664:    Input Parameters:
    1665: .  mat - the matrix

    1667:    Output Parameters:
    1668: +  flag - flag indicating the type of parameters to be returned
    1669:    (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors,
    1670:    MAT_GLOBAL_SUM - sum over all processors)
    1671: -  info - matrix information context

    1673:    Notes:
    1674:    The MatInfo context contains a variety of matrix data, including
    1675:    number of nonzeros allocated and used, number of mallocs during
    1676:    matrix assembly, etc.  Additional information for factored matrices
    1677:    is provided (such as the fill ratio, number of mallocs during
    1678:    factorization, etc.).  Much of this info is printed to STDOUT
    1679:    when using the runtime options
    1680: $       -log_info -mat_view_info

    1682:    Example for C/C++ Users:
    1683:    See the file ${PETSC_DIR}/include/petscmat.h for a complete list of
    1684:    data within the MatInfo context.  For example,
    1685: .vb
    1686:       MatInfo info;
    1687:       Mat     A;
    1688:       double  mal, nz_a, nz_u;

    1690:       MatGetInfo(A,MAT_LOCAL,&info);
    1691:       mal  = info.mallocs;
    1692:       nz_a = info.nz_allocated;
    1693: .ve

    1695:    Example for Fortran Users:
    1696:    Fortran users should declare info as a double precision
    1697:    array of dimension MAT_INFO_SIZE, and then extract the parameters
    1698:    of interest.  See the file ${PETSC_DIR}/include/finclude/petscmat.h
    1699:    a complete list of parameter names.
    1700: .vb
    1701:       double  precision info(MAT_INFO_SIZE)
    1702:       double  precision mal, nz_a
    1703:       Mat     A
    1704:       integer ierr

    1706:       call MatGetInfo(A,MAT_LOCAL,info,ierr)
    1707:       mal = info(MAT_INFO_MALLOCS)
    1708:       nz_a = info(MAT_INFO_NZ_ALLOCATED)
    1709: .ve

    1711:     Level: intermediate

    1713:     Concepts: matrices^getting information on
    1714:
    1715: @*/
    1716: PetscErrorCode PETSCMAT_DLLEXPORT MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info)
    1717: {

    1724:   if (!mat->ops->getinfo) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    1725:   MatPreallocated(mat);
    1726:   (*mat->ops->getinfo)(mat,flag,info);
    1727:   return(0);
    1728: }

    1730: /* ----------------------------------------------------------*/
    1733: /*@C
    1734:    MatILUDTFactor - Performs a drop tolerance ILU factorization.

    1736:    Collective on Mat

    1738:    Input Parameters:
    1739: +  mat - the matrix
    1740: .  row - row permutation
    1741: .  col - column permutation
    1742: -  info - information about the factorization to be done

    1744:    Output Parameters:
    1745: .  fact - the factored matrix

    1747:    Level: developer

    1749:    Notes:
    1750:    Most users should employ the simplified KSP interface for linear solvers
    1751:    instead of working directly with matrix algebra routines such as this.
    1752:    See, e.g., KSPCreate().

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

    1759:     Concepts: matrices^ILUDT factorization

    1761: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
    1762: @*/
    1763: PetscErrorCode PETSCMAT_DLLEXPORT MatILUDTFactor(Mat mat,IS row,IS col,MatFactorInfo *info,Mat *fact)
    1764: {

    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->iludtfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    1777:   MatPreallocated(mat);
    1778:   PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
    1779:   (*mat->ops->iludtfactor)(mat,row,col,info,fact);
    1780:   PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);
    1781:   PetscObjectStateIncrease((PetscObject)*fact);

    1783:   return(0);
    1784: }

    1788: /*@
    1789:    MatLUFactor - Performs in-place LU factorization of matrix.

    1791:    Collective on Mat

    1793:    Input Parameters:
    1794: +  mat - the matrix
    1795: .  row - row permutation
    1796: .  col - column permutation
    1797: -  info - options for factorization, includes
    1798: $          fill - expected fill as ratio of original fill.
    1799: $          dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
    1800: $                   Run with the option -log_info to determine an optimal value to use

    1802:    Notes:
    1803:    Most users should employ the simplified KSP interface for linear solvers
    1804:    instead of working directly with matrix algebra routines such as this.
    1805:    See, e.g., KSPCreate().

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

    1810:    Level: developer

    1812:    Concepts: matrices^LU factorization

    1814: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(),
    1815:           MatGetOrdering(), MatSetUnfactored(), MatFactorInfo

    1817: @*/
    1818: PetscErrorCode PETSCMAT_DLLEXPORT MatLUFactor(Mat mat,IS row,IS col,MatFactorInfo *info)
    1819: {

    1828:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    1829:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1830:   if (!mat->ops->lufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    1831:   MatPreallocated(mat);

    1833:   PetscLogEventBegin(MAT_LUFactor,mat,row,col,0);
    1834:   (*mat->ops->lufactor)(mat,row,col,info);
    1835:   PetscLogEventEnd(MAT_LUFactor,mat,row,col,0);
    1836:   PetscObjectStateIncrease((PetscObject)mat);
    1837:   return(0);
    1838: }

    1842: /*@
    1843:    MatILUFactor - Performs in-place ILU factorization of matrix.

    1845:    Collective on Mat

    1847:    Input Parameters:
    1848: +  mat - the matrix
    1849: .  row - row permutation
    1850: .  col - column permutation
    1851: -  info - structure containing
    1852: $      levels - number of levels of fill.
    1853: $      expected fill - as ratio of original fill.
    1854: $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
    1855:                 missing diagonal entries)

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

    1861:    Most users should employ the simplified KSP interface for linear solvers
    1862:    instead of working directly with matrix algebra routines such as this.
    1863:    See, e.g., KSPCreate().

    1865:    Level: developer

    1867:    Concepts: matrices^ILU factorization

    1869: .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
    1870: @*/
    1871: PetscErrorCode PETSCMAT_DLLEXPORT MatILUFactor(Mat mat,IS row,IS col,MatFactorInfo *info)
    1872: {

    1881:   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
    1882:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    1883:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    1884:   if (!mat->ops->ilufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    1885:   MatPreallocated(mat);

    1887:   PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);
    1888:   (*mat->ops->ilufactor)(mat,row,col,info);
    1889:   PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);
    1890:   PetscObjectStateIncrease((PetscObject)mat);
    1891:   return(0);
    1892: }

    1896: /*@
    1897:    MatLUFactorSymbolic - Performs symbolic LU factorization of matrix.
    1898:    Call this routine before calling MatLUFactorNumeric().

    1900:    Collective on Mat

    1902:    Input Parameters:
    1903: +  mat - the matrix
    1904: .  row, col - row and column permutations
    1905: -  info - options for factorization, includes
    1906: $          fill - expected fill as ratio of original fill.
    1907: $          dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
    1908: $                   Run with the option -log_info to determine an optimal value to use

    1910:    Output Parameter:
    1911: .  fact - new matrix that has been symbolically factored

    1913:    Notes:
    1914:    See the users manual for additional information about
    1915:    choosing the fill factor for better efficiency.

    1917:    Most users should employ the simplified KSP interface for linear solvers
    1918:    instead of working directly with matrix algebra routines such as this.
    1919:    See, e.g., KSPCreate().

    1921:    Level: developer

    1923:    Concepts: matrices^LU symbolic factorization

    1925: .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
    1926: @*/
    1927: PetscErrorCode PETSCMAT_DLLEXPORT MatLUFactorSymbolic(Mat mat,IS row,IS col,MatFactorInfo *info,Mat *fact)
    1928: {

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

    1943:   PetscLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);
    1944:   (*mat->ops->lufactorsymbolic)(mat,row,col,info,fact);
    1945:   PetscLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);
    1946:   PetscObjectStateIncrease((PetscObject)*fact);
    1947:   return(0);
    1948: }

    1952: /*@
    1953:    MatLUFactorNumeric - Performs numeric LU factorization of a matrix.
    1954:    Call this routine after first calling MatLUFactorSymbolic().

    1956:    Collective on Mat

    1958:    Input Parameters:
    1959: +  mat - the matrix
    1960: .  info - options for factorization
    1961: -  fact - the matrix generated for the factor, from MatLUFactorSymbolic()

    1963:    Notes:
    1964:    See MatLUFactor() for in-place factorization.  See
    1965:    MatCholeskyFactorNumeric() for the symmetric, positive definite case.

    1967:    Most users should employ the simplified KSP interface for linear solvers
    1968:    instead of working directly with matrix algebra routines such as this.
    1969:    See, e.g., KSPCreate().

    1971:    Level: developer

    1973:    Concepts: matrices^LU numeric factorization

    1975: .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor()
    1976: @*/
    1977: PetscErrorCode PETSCMAT_DLLEXPORT MatLUFactorNumeric(Mat mat,MatFactorInfo *info,Mat *fact)
    1978: {

    1986:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    1987:   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
    1988:     SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dimensions are different %D should = %D %D should = %D",
    1989:             mat->M,(*fact)->M,mat->N,(*fact)->N);
    1990:   }
    1991:   if (!(*fact)->ops->lufactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    1992:   MatPreallocated(mat);
    1993:   PetscLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0);
    1994:   (*(*fact)->ops->lufactornumeric)(mat,info,fact);
    1995:   PetscLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0);

    1997:   MatView_Private(*fact);
    1998:   PetscObjectStateIncrease((PetscObject)*fact);
    1999:   return(0);
    2000: }

    2004: /*@
    2005:    MatCholeskyFactor - Performs in-place Cholesky factorization of a
    2006:    symmetric matrix.

    2008:    Collective on Mat

    2010:    Input Parameters:
    2011: +  mat - the matrix
    2012: .  perm - row and column permutations
    2013: -  f - expected fill as ratio of original fill

    2015:    Notes:
    2016:    See MatLUFactor() for the nonsymmetric case.  See also
    2017:    MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric().

    2019:    Most users should employ the simplified KSP interface for linear solvers
    2020:    instead of working directly with matrix algebra routines such as this.
    2021:    See, e.g., KSPCreate().

    2023:    Level: developer

    2025:    Concepts: matrices^Cholesky factorization

    2027: .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric()
    2028:           MatGetOrdering()

    2030: @*/
    2031: PetscErrorCode PETSCMAT_DLLEXPORT MatCholeskyFactor(Mat mat,IS perm,MatFactorInfo *info)
    2032: {

    2040:   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
    2041:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    2042:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2043:   if (!mat->ops->choleskyfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    2044:   MatPreallocated(mat);

    2046:   PetscLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);
    2047:   (*mat->ops->choleskyfactor)(mat,perm,info);
    2048:   PetscLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);
    2049:   PetscObjectStateIncrease((PetscObject)mat);
    2050:   return(0);
    2051: }

    2055: /*@
    2056:    MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization
    2057:    of a symmetric matrix.

    2059:    Collective on Mat

    2061:    Input Parameters:
    2062: +  mat - the matrix
    2063: .  perm - row and column permutations
    2064: -  info - options for factorization, includes
    2065: $          fill - expected fill as ratio of original fill.
    2066: $          dtcol - pivot tolerance (0 no pivot, 1 full column pivoting)
    2067: $                   Run with the option -log_info to determine an optimal value to use

    2069:    Output Parameter:
    2070: .  fact - the factored matrix

    2072:    Notes:
    2073:    See MatLUFactorSymbolic() for the nonsymmetric case.  See also
    2074:    MatCholeskyFactor() and MatCholeskyFactorNumeric().

    2076:    Most users should employ the simplified KSP interface for linear solvers
    2077:    instead of working directly with matrix algebra routines such as this.
    2078:    See, e.g., KSPCreate().

    2080:    Level: developer

    2082:    Concepts: matrices^Cholesky symbolic factorization

    2084: .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric()
    2085:           MatGetOrdering()

    2087: @*/
    2088: PetscErrorCode PETSCMAT_DLLEXPORT MatCholeskyFactorSymbolic(Mat mat,IS perm,MatFactorInfo *info,Mat *fact)
    2089: {

    2098:   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square");
    2099:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    2100:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2101:   if (!mat->ops->choleskyfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    2102:   MatPreallocated(mat);

    2104:   PetscLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
    2105:   (*mat->ops->choleskyfactorsymbolic)(mat,perm,info,fact);
    2106:   PetscLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);
    2107:   PetscObjectStateIncrease((PetscObject)*fact);
    2108:   return(0);
    2109: }

    2113: /*@
    2114:    MatCholeskyFactorNumeric - Performs numeric Cholesky factorization
    2115:    of a symmetric matrix. Call this routine after first calling
    2116:    MatCholeskyFactorSymbolic().

    2118:    Collective on Mat

    2120:    Input Parameter:
    2121: .  mat - the initial matrix
    2122: .  info - options for factorization
    2123: -  fact - the symbolic factor of mat

    2125:    Output Parameter:
    2126: .  fact - the factored matrix

    2128:    Notes:
    2129:    Most users should employ the simplified KSP interface for linear solvers
    2130:    instead of working directly with matrix algebra routines such as this.
    2131:    See, e.g., KSPCreate().

    2133:    Level: developer

    2135:    Concepts: matrices^Cholesky numeric factorization

    2137: .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric()
    2138: @*/
    2139: PetscErrorCode PETSCMAT_DLLEXPORT MatCholeskyFactorNumeric(Mat mat,MatFactorInfo *info,Mat *fact)
    2140: {

    2147:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    2148:   if (!(*fact)->ops->choleskyfactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    2149:   if (mat->M != (*fact)->M || mat->N != (*fact)->N) {
    2150:     SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dim %D should = %D %D should = %D",mat->M,(*fact)->M,mat->N,(*fact)->N);
    2151:   }
    2152:   MatPreallocated(mat);

    2154:   PetscLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
    2155:   (*(*fact)->ops->choleskyfactornumeric)(mat,info,fact);
    2156:   PetscLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0);
    2157:   PetscObjectStateIncrease((PetscObject)*fact);
    2158:   return(0);
    2159: }

    2161: /* ----------------------------------------------------------------*/
    2164: /*@
    2165:    MatSolve - Solves A x = b, given a factored matrix.

    2167:    Collective on Mat and Vec

    2169:    Input Parameters:
    2170: +  mat - the factored matrix
    2171: -  b - the right-hand-side vector

    2173:    Output Parameter:
    2174: .  x - the result vector

    2176:    Notes:
    2177:    The vectors b and x cannot be the same.  I.e., one cannot
    2178:    call MatSolve(A,x,x).

    2180:    Notes:
    2181:    Most users should employ the simplified KSP interface for linear solvers
    2182:    instead of working directly with matrix algebra routines such as this.
    2183:    See, e.g., KSPCreate().

    2185:    Level: developer

    2187:    Concepts: matrices^triangular solves

    2189: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd()
    2190: @*/
    2191: PetscErrorCode PETSCMAT_DLLEXPORT MatSolve(Mat mat,Vec b,Vec x)
    2192: {

    2202:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
    2203:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
    2204:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->N,x->N);
    2205:   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->M,b->N);
    2206:   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->m,b->n);
    2207:   if (!mat->M && !mat->N) return(0);
    2208:   if (!mat->ops->solve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    2209:   MatPreallocated(mat);

    2211:   PetscLogEventBegin(MAT_Solve,mat,b,x,0);
    2212:   (*mat->ops->solve)(mat,b,x);
    2213:   PetscLogEventEnd(MAT_Solve,mat,b,x,0);
    2214:   PetscObjectStateIncrease((PetscObject)x);
    2215:   return(0);
    2216: }

    2220: /* @
    2221:    MatForwardSolve - Solves L x = b, given a factored matrix, A = LU.

    2223:    Collective on Mat and Vec

    2225:    Input Parameters:
    2226: +  mat - the factored matrix
    2227: -  b - the right-hand-side vector

    2229:    Output Parameter:
    2230: .  x - the result vector

    2232:    Notes:
    2233:    MatSolve() should be used for most applications, as it performs
    2234:    a forward solve followed by a backward solve.

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

    2239:    Most users should employ the simplified KSP interface for linear solvers
    2240:    instead of working directly with matrix algebra routines such as this.
    2241:    See, e.g., KSPCreate().

    2243:    Level: developer

    2245:    Concepts: matrices^forward solves

    2247: .seealso: MatSolve(), MatBackwardSolve()
    2248: @ */
    2249: PetscErrorCode PETSCMAT_DLLEXPORT MatForwardSolve(Mat mat,Vec b,Vec x)
    2250: {

    2260:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
    2261:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
    2262:   if (!mat->ops->forwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    2263:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->N,x->N);
    2264:   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->M,b->N);
    2265:   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->m,b->n);
    2266:   MatPreallocated(mat);
    2267:   PetscLogEventBegin(MAT_ForwardSolve,mat,b,x,0);
    2268:   (*mat->ops->forwardsolve)(mat,b,x);
    2269:   PetscLogEventEnd(MAT_ForwardSolve,mat,b,x,0);
    2270:   PetscObjectStateIncrease((PetscObject)x);
    2271:   return(0);
    2272: }

    2276: /* @
    2277:    MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU.

    2279:    Collective on Mat and Vec

    2281:    Input Parameters:
    2282: +  mat - the factored matrix
    2283: -  b - the right-hand-side vector

    2285:    Output Parameter:
    2286: .  x - the result vector

    2288:    Notes:
    2289:    MatSolve() should be used for most applications, as it performs
    2290:    a forward solve followed by a backward solve.

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

    2295:    Most users should employ the simplified KSP interface for linear solvers
    2296:    instead of working directly with matrix algebra routines such as this.
    2297:    See, e.g., KSPCreate().

    2299:    Level: developer

    2301:    Concepts: matrices^backward solves

    2303: .seealso: MatSolve(), MatForwardSolve()
    2304: @ */
    2305: PetscErrorCode PETSCMAT_DLLEXPORT MatBackwardSolve(Mat mat,Vec b,Vec x)
    2306: {

    2316:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
    2317:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
    2318:   if (!mat->ops->backwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    2319:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->N,x->N);
    2320:   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->M,b->N);
    2321:   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->m,b->n);
    2322:   MatPreallocated(mat);

    2324:   PetscLogEventBegin(MAT_BackwardSolve,mat,b,x,0);
    2325:   (*mat->ops->backwardsolve)(mat,b,x);
    2326:   PetscLogEventEnd(MAT_BackwardSolve,mat,b,x,0);
    2327:   PetscObjectStateIncrease((PetscObject)x);
    2328:   return(0);
    2329: }

    2333: /*@
    2334:    MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix.

    2336:    Collective on Mat and Vec

    2338:    Input Parameters:
    2339: +  mat - the factored matrix
    2340: .  b - the right-hand-side vector
    2341: -  y - the vector to be added to

    2343:    Output Parameter:
    2344: .  x - the result vector

    2346:    Notes:
    2347:    The vectors b and x cannot be the same.  I.e., one cannot
    2348:    call MatSolveAdd(A,x,y,x).

    2350:    Most users should employ the simplified KSP interface for linear solvers
    2351:    instead of working directly with matrix algebra routines such as this.
    2352:    See, e.g., KSPCreate().

    2354:    Level: developer

    2356:    Concepts: matrices^triangular solves

    2358: .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd()
    2359: @*/
    2360: PetscErrorCode PETSCMAT_DLLEXPORT MatSolveAdd(Mat mat,Vec b,Vec y,Vec x)
    2361: {
    2362:   PetscScalar    one = 1.0;
    2363:   Vec            tmp;

    2375:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
    2376:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
    2377:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->N,x->N);
    2378:   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->M,b->N);
    2379:   if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->M,y->N);
    2380:   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->m,b->n);
    2381:   if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %D %D",x->n,y->n);
    2382:   MatPreallocated(mat);

    2384:   PetscLogEventBegin(MAT_SolveAdd,mat,b,x,y);
    2385:   if (mat->ops->solveadd)  {
    2386:     (*mat->ops->solveadd)(mat,b,y,x);
    2387:   } else {
    2388:     /* do the solve then the add manually */
    2389:     if (x != y) {
    2390:       MatSolve(mat,b,x);
    2391:       VecAXPY(x,one,y);
    2392:     } else {
    2393:       VecDuplicate(x,&tmp);
    2394:       PetscLogObjectParent(mat,tmp);
    2395:       VecCopy(x,tmp);
    2396:       MatSolve(mat,b,x);
    2397:       VecAXPY(x,one,tmp);
    2398:       VecDestroy(tmp);
    2399:     }
    2400:   }
    2401:   PetscLogEventEnd(MAT_SolveAdd,mat,b,x,y);
    2402:   PetscObjectStateIncrease((PetscObject)x);
    2403:   return(0);
    2404: }

    2408: /*@
    2409:    MatSolveTranspose - Solves A' x = b, given a factored matrix.

    2411:    Collective on Mat and Vec

    2413:    Input Parameters:
    2414: +  mat - the factored matrix
    2415: -  b - the right-hand-side vector

    2417:    Output Parameter:
    2418: .  x - the result vector

    2420:    Notes:
    2421:    The vectors b and x cannot be the same.  I.e., one cannot
    2422:    call MatSolveTranspose(A,x,x).

    2424:    Most users should employ the simplified KSP interface for linear solvers
    2425:    instead of working directly with matrix algebra routines such as this.
    2426:    See, e.g., KSPCreate().

    2428:    Level: developer

    2430:    Concepts: matrices^triangular solves

    2432: .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd()
    2433: @*/
    2434: PetscErrorCode PETSCMAT_DLLEXPORT MatSolveTranspose(Mat mat,Vec b,Vec x)
    2435: {

    2445:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
    2446:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
    2447:   if (!mat->ops->solvetranspose) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s",mat->type_name);
    2448:   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->M,x->N);
    2449:   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->N,b->N);
    2450:   MatPreallocated(mat);
    2451:   PetscLogEventBegin(MAT_SolveTranspose,mat,b,x,0);
    2452:   (*mat->ops->solvetranspose)(mat,b,x);
    2453:   PetscLogEventEnd(MAT_SolveTranspose,mat,b,x,0);
    2454:   PetscObjectStateIncrease((PetscObject)x);
    2455:   return(0);
    2456: }

    2460: /*@
    2461:    MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a
    2462:                       factored matrix.

    2464:    Collective on Mat and Vec

    2466:    Input Parameters:
    2467: +  mat - the factored matrix
    2468: .  b - the right-hand-side vector
    2469: -  y - the vector to be added to

    2471:    Output Parameter:
    2472: .  x - the result vector

    2474:    Notes:
    2475:    The vectors b and x cannot be the same.  I.e., one cannot
    2476:    call MatSolveTransposeAdd(A,x,y,x).

    2478:    Most users should employ the simplified KSP interface for linear solvers
    2479:    instead of working directly with matrix algebra routines such as this.
    2480:    See, e.g., KSPCreate().

    2482:    Level: developer

    2484:    Concepts: matrices^triangular solves

    2486: .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose()
    2487: @*/
    2488: PetscErrorCode PETSCMAT_DLLEXPORT MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x)
    2489: {
    2490:   PetscScalar    one = 1.0;
    2492:   Vec            tmp;

    2503:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
    2504:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
    2505:   if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->M,x->N);
    2506:   if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->N,b->N);
    2507:   if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %D %D",mat->N,y->N);
    2508:   if (x->n != y->n)   SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %D %D",x->n,y->n);
    2509:   MatPreallocated(mat);

    2511:   PetscLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);
    2512:   if (mat->ops->solvetransposeadd) {
    2513:     (*mat->ops->solvetransposeadd)(mat,b,y,x);
    2514:   } else {
    2515:     /* do the solve then the add manually */
    2516:     if (x != y) {
    2517:       MatSolveTranspose(mat,b,x);
    2518:       VecAXPY(x,one,y);
    2519:     } else {
    2520:       VecDuplicate(x,&tmp);
    2521:       PetscLogObjectParent(mat,tmp);
    2522:       VecCopy(x,tmp);
    2523:       MatSolveTranspose(mat,b,x);
    2524:       VecAXPY(x,one,tmp);
    2525:       VecDestroy(tmp);
    2526:     }
    2527:   }
    2528:   PetscLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);
    2529:   PetscObjectStateIncrease((PetscObject)x);
    2530:   return(0);
    2531: }
    2532: /* ----------------------------------------------------------------*/

    2536: /*@
    2537:    MatRelax - Computes relaxation (SOR, Gauss-Seidel) sweeps.

    2539:    Collective on Mat and Vec

    2541:    Input Parameters:
    2542: +  mat - the matrix
    2543: .  b - the right hand side
    2544: .  omega - the relaxation factor
    2545: .  flag - flag indicating the type of SOR (see below)
    2546: .  shift -  diagonal shift
    2547: -  its - the number of iterations
    2548: -  lits - the number of local iterations

    2550:    Output Parameters:
    2551: .  x - the solution (can contain an initial guess)

    2553:    SOR Flags:
    2554: .     SOR_FORWARD_SWEEP - forward SOR
    2555: .     SOR_BACKWARD_SWEEP - backward SOR
    2556: .     SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR)
    2557: .     SOR_LOCAL_FORWARD_SWEEP - local forward SOR
    2558: .     SOR_LOCAL_BACKWARD_SWEEP - local forward SOR
    2559: .     SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR
    2560: .     SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies
    2561:          upper/lower triangular part of matrix to
    2562:          vector (with omega)
    2563: .     SOR_ZERO_INITIAL_GUESS - zero initial guess

    2565:    Notes:
    2566:    SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and
    2567:    SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings
    2568:    on each processor.

    2570:    Application programmers will not generally use MatRelax() directly,
    2571:    but instead will employ the KSP/PC interface.

    2573:    Notes for Advanced Users:
    2574:    The flags are implemented as bitwise inclusive or operations.
    2575:    For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP)
    2576:    to specify a zero initial guess for SSOR.

    2578:    Most users should employ the simplified KSP interface for linear solvers
    2579:    instead of working directly with matrix algebra routines such as this.
    2580:    See, e.g., KSPCreate().

    2582:    See also, MatPBRelax(). This routine will automatically call the point block
    2583:    version if the point version is not available.

    2585:    Level: developer

    2587:    Concepts: matrices^relaxation
    2588:    Concepts: matrices^SOR
    2589:    Concepts: matrices^Gauss-Seidel

    2591: @*/
    2592: PetscErrorCode PETSCMAT_DLLEXPORT MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec x)
    2593: {

    2603:   if (!mat->ops->relax && !mat->ops->pbrelax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    2604:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    2605:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2606:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->N,x->N);
    2607:   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->M,b->N);
    2608:   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->m,b->n);
    2609:   MatPreallocated(mat);
    2610:   PetscLogEventBegin(MAT_Relax,mat,b,x,0);
    2611:   if (mat->ops->relax) {
    2612:     ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,lits,x);
    2613:   } else {
    2614:     ierr =(*mat->ops->pbrelax)(mat,b,omega,flag,shift,its,lits,x);
    2615:   }
    2616:   PetscLogEventEnd(MAT_Relax,mat,b,x,0);
    2617:   PetscObjectStateIncrease((PetscObject)x);
    2618:   return(0);
    2619: }

    2623: /*@
    2624:    MatPBRelax - Computes relaxation (SOR, Gauss-Seidel) sweeps.

    2626:    Collective on Mat and Vec

    2628:    See MatRelax() for usage

    2630:    For multi-component PDEs where the Jacobian is stored in a point block format
    2631:    (with the PETSc BAIJ matrix formats) the relaxation is done one point block at
    2632:    a time. That is, the small (for example, 4 by 4) blocks along the diagonal are solved
    2633:    simultaneously (that is a 4 by 4 linear solve is done) to update all the values at a point.

    2635:    Level: developer

    2637: @*/
    2638: PetscErrorCode PETSCMAT_DLLEXPORT MatPBRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec x)
    2639: {

    2649:   if (!mat->ops->pbrelax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    2650:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    2651:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2652:   if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %D %D",mat->N,x->N);
    2653:   if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %D %D",mat->M,b->N);
    2654:   if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %D %D",mat->m,b->n);
    2655:   MatPreallocated(mat);

    2657:   PetscLogEventBegin(MAT_Relax,mat,b,x,0);
    2658:   ierr =(*mat->ops->pbrelax)(mat,b,omega,flag,shift,its,lits,x);
    2659:   PetscLogEventEnd(MAT_Relax,mat,b,x,0);
    2660:   PetscObjectStateIncrease((PetscObject)x);
    2661:   return(0);
    2662: }

    2666: /*
    2667:       Default matrix copy routine.
    2668: */
    2669: PetscErrorCode MatCopy_Basic(Mat A,Mat B,MatStructure str)
    2670: {
    2671:   PetscErrorCode    ierr;
    2672:   PetscInt          i,rstart,rend,nz;
    2673:   const PetscInt    *cwork;
    2674:   const PetscScalar *vwork;

    2677:   if (B->assembled) {
    2678:     MatZeroEntries(B);
    2679:   }
    2680:   MatGetOwnershipRange(A,&rstart,&rend);
    2681:   for (i=rstart; i<rend; i++) {
    2682:     MatGetRow(A,i,&nz,&cwork,&vwork);
    2683:     MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);
    2684:     MatRestoreRow(A,i,&nz,&cwork,&vwork);
    2685:   }
    2686:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    2687:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    2688:   PetscObjectStateIncrease((PetscObject)B);
    2689:   return(0);
    2690: }

    2694: /*@C
    2695:    MatCopy - Copys a matrix to another matrix.

    2697:    Collective on Mat

    2699:    Input Parameters:
    2700: +  A - the matrix
    2701: -  str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN

    2703:    Output Parameter:
    2704: .  B - where the copy is put

    2706:    Notes:
    2707:    If you use SAME_NONZERO_PATTERN then the two matrices had better have the
    2708:    same nonzero pattern or the routine will crash.

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

    2714:    Level: intermediate
    2715:
    2716:    Concepts: matrices^copying

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

    2720: @*/
    2721: PetscErrorCode PETSCMAT_DLLEXPORT MatCopy(Mat A,Mat B,MatStructure str)
    2722: {

    2730:   MatPreallocated(B);
    2732:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    2733:   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2734:   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);
    2735:   MatPreallocated(A);

    2737:   PetscLogEventBegin(MAT_Copy,A,B,0,0);
    2738:   if (A->ops->copy) {
    2739:     (*A->ops->copy)(A,B,str);
    2740:   } else { /* generic conversion */
    2741:     MatCopy_Basic(A,B,str);
    2742:   }
    2743:   if (A->mapping) {
    2744:     if (B->mapping) {ISLocalToGlobalMappingDestroy(B->mapping);B->mapping = 0;}
    2745:     MatSetLocalToGlobalMapping(B,A->mapping);
    2746:   }
    2747:   if (A->bmapping) {
    2748:     if (B->bmapping) {ISLocalToGlobalMappingDestroy(B->bmapping);B->bmapping = 0;}
    2749:     MatSetLocalToGlobalMappingBlock(B,A->mapping);
    2750:   }
    2751:   PetscLogEventEnd(MAT_Copy,A,B,0,0);
    2752:   PetscObjectStateIncrease((PetscObject)B);
    2753:   return(0);
    2754: }

    2756:  #include petscsys.h
    2757: PetscTruth MatConvertRegisterAllCalled = PETSC_FALSE;
    2758: PetscFList MatConvertList              = 0;

    2762: /*@C
    2763:     MatConvertRegister - Allows one to register a routine that converts a sparse matrix
    2764:         from one format to another.

    2766:   Not Collective

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

    2772:   Level: developer

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

    2776: @*/
    2777: PetscErrorCode PETSCMAT_DLLEXPORT MatConvertRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat*))
    2778: {
    2780:   char           fullname[PETSC_MAX_PATH_LEN];

    2783:   PetscFListConcat(path,name,fullname);
    2784:   PetscFListAdd(&MatConvertList,sname,fullname,(void (*)(void))function);
    2785:   return(0);
    2786: }

    2790: /*@C
    2791:    MatConvert - Converts a matrix to another matrix, either of the same
    2792:    or different type.

    2794:    Collective on Mat

    2796:    Input Parameters:
    2797: +  mat - the matrix
    2798: .  newtype - new matrix type.  Use MATSAME to create a new matrix of the
    2799:    same type as the original matrix.
    2800: -  reuse - denotes if the destination matrix is to be created or reused.  Currently
    2801:    MAT_REUSE_MATRIX is only supported for inplace conversion, otherwise use
    2802:    MAT_INITIAL_MATRIX.
    2803:    Output Parameter:
    2804: .  M - pointer to place new matrix

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

    2811:    Level: intermediate

    2813:    Concepts: matrices^converting between storage formats

    2815: .seealso: MatCopy(), MatDuplicate()
    2816: @*/
    2817: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert(Mat mat, MatType newtype,MatReuse reuse,Mat *M)
    2818: {
    2819:   PetscErrorCode         ierr;
    2820:   PetscTruth             sametype,issame,flg;
    2821:   char                   convname[256],mtype[256];
    2822:   Mat                    B;

    2828:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    2829:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2830:   MatPreallocated(mat);

    2832:   PetscOptionsGetString(PETSC_NULL,"-matconvert_type",mtype,256,&flg);
    2833:   if (flg) {
    2834:     newtype = mtype;
    2835:   }
    2836:   PetscLogEventBegin(MAT_Convert,mat,0,0,0);
    2837:
    2838:   PetscTypeCompare((PetscObject)mat,newtype,&sametype);
    2839:   PetscStrcmp(newtype,"same",&issame);
    2840:   if ((reuse==MAT_REUSE_MATRIX) && (mat != *M)) {
    2841:     SETERRQ(PETSC_ERR_SUP,"MAT_REUSE_MATRIX only supported for inplace convertion currently");
    2842:   }
    2843:   if ((sametype || issame) && (reuse==MAT_INITIAL_MATRIX) && mat->ops->duplicate) {
    2844:     (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);
    2845:   } else {
    2846:     PetscErrorCode (*conv)(Mat, MatType,MatReuse,Mat*)=PETSC_NULL;
    2847:     /*
    2848:        Order of precedence:
    2849:        1) See if a specialized converter is known to the current matrix.
    2850:        2) See if a specialized converter is known to the desired matrix class.
    2851:        3) See if a good general converter is registered for the desired class
    2852:           (as of 6/27/03 only MATMPIADJ falls into this category).
    2853:        4) See if a good general converter is known for the current matrix.
    2854:        5) Use a really basic converter.
    2855:     */
    2856:     PetscStrcpy(convname,"MatConvert_");
    2857:     PetscStrcat(convname,mat->type_name);
    2858:     PetscStrcat(convname,"_");
    2859:     PetscStrcat(convname,newtype);
    2860:     PetscStrcat(convname,"_C");
    2861:     PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);

    2863:     if (!conv) {
    2864:       MatCreate(mat->comm,&B);
    2865:       MatSetSizes(B,0,0,0,0);
    2866:       MatSetType(B,newtype);
    2867:       PetscObjectQueryFunction((PetscObject)B,convname,(void (**)(void))&conv);
    2868:       MatDestroy(B);
    2869:       if (!conv) {
    2870:         if (!MatConvertRegisterAllCalled) {
    2871:           MatConvertRegisterAll(PETSC_NULL);
    2872:         }
    2873:         PetscFListFind(mat->comm,MatConvertList,newtype,(void(**)(void))&conv);
    2874:         if (!conv) {
    2875:           if (mat->ops->convert) {
    2876:             conv = mat->ops->convert;
    2877:           } else {
    2878:             conv = MatConvert_Basic;
    2879:           }
    2880:         }
    2881:       }
    2882:     }
    2883:     (*conv)(mat,newtype,reuse,M);
    2884:   }
    2885:   B = *M;
    2886:   PetscLogEventEnd(MAT_Convert,mat,0,0,0);
    2887:   PetscObjectStateIncrease((PetscObject)B);
    2888:   return(0);
    2889: }

    2894: /*@C
    2895:    MatDuplicate - Duplicates a matrix including the non-zero structure.

    2897:    Collective on Mat

    2899:    Input Parameters:
    2900: +  mat - the matrix
    2901: -  op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero
    2902:         values as well or not

    2904:    Output Parameter:
    2905: .  M - pointer to place new matrix

    2907:    Level: intermediate

    2909:    Concepts: matrices^duplicating

    2911: .seealso: MatCopy(), MatConvert()
    2912: @*/
    2913: PetscErrorCode PETSCMAT_DLLEXPORT MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
    2914: {
    2916:   Mat            B;

    2922:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    2923:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    2924:   MatPreallocated(mat);

    2926:   *M  = 0;
    2927:   PetscLogEventBegin(MAT_Convert,mat,0,0,0);
    2928:   if (!mat->ops->duplicate) {
    2929:     SETERRQ(PETSC_ERR_SUP,"Not written for this matrix type");
    2930:   }
    2931:   (*mat->ops->duplicate)(mat,op,M);
    2932:   B = *M;
    2933:   if (mat->mapping) {
    2934:     MatSetLocalToGlobalMapping(B,mat->mapping);
    2935:   }
    2936:   if (mat->bmapping) {
    2937:     MatSetLocalToGlobalMappingBlock(B,mat->bmapping);
    2938:   }
    2939:   if (mat->rmap){
    2940:     if (!B->rmap){
    2941:       PetscMapCreateMPI(B->comm,B->m,B->M,&B->rmap);
    2942:     }
    2943:     PetscMemcpy(B->rmap,mat->rmap,sizeof(PetscMap));
    2944:   }
    2945:   if (mat->cmap){
    2946:     if (!B->cmap){
    2947:       PetscMapCreateMPI(B->comm,B->n,B->N,&B->cmap);
    2948:     }
    2949:     PetscMemcpy(B->cmap,mat->cmap,sizeof(PetscMap));
    2950:   }
    2951:   PetscLogEventEnd(MAT_Convert,mat,0,0,0);
    2952:   PetscObjectStateIncrease((PetscObject)B);
    2953:   return(0);
    2954: }

    2958: /*@
    2959:    MatGetDiagonal - Gets the diagonal of a matrix.

    2961:    Collective on Mat and Vec

    2963:    Input Parameters:
    2964: +  mat - the matrix
    2965: -  v - the vector for storing the diagonal

    2967:    Output Parameter:
    2968: .  v - the diagonal of the matrix

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

    2976:    Level: intermediate

    2978:    Concepts: matrices^accessing diagonals

    2980: .seealso: MatGetRow(), MatGetSubmatrices(), MatGetSubmatrix(), MatGetRowMax()
    2981: @*/
    2982: PetscErrorCode PETSCMAT_DLLEXPORT MatGetDiagonal(Mat mat,Vec v)
    2983: {

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

    2994:   (*mat->ops->getdiagonal)(mat,v);
    2995:   PetscObjectStateIncrease((PetscObject)v);
    2996:   return(0);
    2997: }

    3001: /*@
    3002:    MatGetRowMax - Gets the maximum value (in absolute value) of each
    3003:         row of the matrix

    3005:    Collective on Mat and Vec

    3007:    Input Parameters:
    3008: .  mat - the matrix

    3010:    Output Parameter:
    3011: .  v - the vector for storing the maximums

    3013:    Level: intermediate

    3015:    Concepts: matrices^getting row maximums

    3017: .seealso: MatGetDiagonal(), MatGetSubmatrices(), MatGetSubmatrix()
    3018: @*/
    3019: PetscErrorCode PETSCMAT_DLLEXPORT MatGetRowMax(Mat mat,Vec v)
    3020: {

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

    3031:   (*mat->ops->getrowmax)(mat,v);
    3032:   PetscObjectStateIncrease((PetscObject)v);
    3033:   return(0);
    3034: }

    3038: /*@C
    3039:    MatTranspose - Computes an in-place or out-of-place transpose of a matrix.

    3041:    Collective on Mat

    3043:    Input Parameter:
    3044: .  mat - the matrix to transpose

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

    3049:    Level: intermediate

    3051:    Concepts: matrices^transposing

    3053: .seealso: MatMultTranspose(), MatMultTransposeAdd(), MatIsTranspose()
    3054: @*/
    3055: PetscErrorCode PETSCMAT_DLLEXPORT MatTranspose(Mat mat,Mat *B)
    3056: {

    3062:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    3063:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    3064:   if (!mat->ops->transpose) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    3065:   MatPreallocated(mat);

    3067:   PetscLogEventBegin(MAT_Transpose,mat,0,0,0);
    3068:   (*mat->ops->transpose)(mat,B);
    3069:   PetscLogEventEnd(MAT_Transpose,mat,0,0,0);
    3070:   if (B) {PetscObjectStateIncrease((PetscObject)*B);}
    3071:   return(0);
    3072: }

    3076: /*@C
    3077:    MatIsTranspose - Test whether a matrix is another one's transpose,
    3078:         or its own, in which case it tests symmetry.

    3080:    Collective on Mat

    3082:    Input Parameter:
    3083: +  A - the matrix to test
    3084: -  B - the matrix to test against, this can equal the first parameter

    3086:    Output Parameters:
    3087: .  flg - the result

    3089:    Notes:
    3090:    Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm
    3091:    has a running time of the order of the number of nonzeros; the parallel
    3092:    test involves parallel copies of the block-offdiagonal parts of the matrix.

    3094:    Level: intermediate

    3096:    Concepts: matrices^transposing, matrix^symmetry

    3098: .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian()
    3099: @*/
    3100: PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose(Mat A,Mat B,PetscReal tol,PetscTruth *flg)
    3101: {
    3102:   PetscErrorCode ierr,(*f)(Mat,Mat,PetscReal,PetscTruth*),(*g)(Mat,Mat,PetscReal,PetscTruth*);

    3108:   PetscObjectQueryFunction((PetscObject)A,"MatIsTranspose_C",(void (**)(void))&f);
    3109:   PetscObjectQueryFunction((PetscObject)B,"MatIsTranspose_C",(void (**)(void))&g);
    3110:   if (f && g) {
    3111:     if (f==g) {
    3112:       (*f)(A,B,tol,flg);
    3113:     } else {
    3114:       SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Matrices do not have the same comparator for symmetry test");
    3115:     }
    3116:   }
    3117:   return(0);
    3118: }

    3122: /*@C
    3123:    MatPermute - Creates a new matrix with rows and columns permuted from the
    3124:    original.

    3126:    Collective on Mat

    3128:    Input Parameters:
    3129: +  mat - the matrix to permute
    3130: .  row - row permutation, each processor supplies only the permutation for its rows
    3131: -  col - column permutation, each processor needs the entire column permutation, that is
    3132:          this is the same size as the total number of columns in the matrix

    3134:    Output Parameters:
    3135: .  B - the permuted matrix

    3137:    Level: advanced

    3139:    Concepts: matrices^permuting

    3141: .seealso: MatGetOrdering()
    3142: @*/
    3143: PetscErrorCode PETSCMAT_DLLEXPORT MatPermute(Mat mat,IS row,IS col,Mat *B)
    3144: {

    3153:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    3154:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    3155:   if (!mat->ops->permute) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    3156:   MatPreallocated(mat);

    3158:   (*mat->ops->permute)(mat,row,col,B);
    3159:   PetscObjectStateIncrease((PetscObject)*B);
    3160:   return(0);
    3161: }

    3165: /*@C
    3166:   MatPermuteSparsify - Creates a new matrix with rows and columns permuted from the
    3167:   original and sparsified to the prescribed tolerance.

    3169:   Collective on Mat

    3171:   Input Parameters:
    3172: + A    - The matrix to permute
    3173: . band - The half-bandwidth of the sparsified matrix, or PETSC_DECIDE
    3174: . frac - The half-bandwidth as a fraction of the total size, or 0.0
    3175: . tol  - The drop tolerance
    3176: . rowp - The row permutation
    3177: - colp - The column permutation

    3179:   Output Parameter:
    3180: . B    - The permuted, sparsified matrix

    3182:   Level: advanced

    3184:   Note:
    3185:   The default behavior (band = PETSC_DECIDE and frac = 0.0) is to
    3186:   restrict the half-bandwidth of the resulting matrix to 5% of the
    3187:   total matrix size.

    3189: .keywords: matrix, permute, sparsify

    3191: .seealso: MatGetOrdering(), MatPermute()
    3192: @*/
    3193: PetscErrorCode PETSCMAT_DLLEXPORT MatPermuteSparsify(Mat A, PetscInt band, PetscReal frac, PetscReal tol, IS rowp, IS colp, Mat *B)
    3194: {
    3195:   IS                irowp, icolp;
    3196:   PetscInt          *rows, *cols;
    3197:   PetscInt          M, N, locRowStart, locRowEnd;
    3198:   PetscInt          nz, newNz;
    3199:   const PetscInt    *cwork;
    3200:   PetscInt          *cnew;
    3201:   const PetscScalar *vwork;
    3202:   PetscScalar       *vnew;
    3203:   PetscInt          bw, issize;
    3204:   PetscInt          row, locRow, newRow, col, newCol;
    3205:   PetscErrorCode    ierr;

    3212:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
    3213:   if (A->factor)     SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
    3214:   if (!A->ops->permutesparsify) {
    3215:     MatGetSize(A, &M, &N);
    3216:     MatGetOwnershipRange(A, &locRowStart, &locRowEnd);
    3217:     ISGetSize(rowp, &issize);
    3218:     if (issize != M) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %D for row permutation, should be %D", issize, M);
    3219:     ISGetSize(colp, &issize);
    3220:     if (issize != N) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %D for column permutation, should be %D", issize, N);
    3221:     ISInvertPermutation(rowp, 0, &irowp);
    3222:     ISGetIndices(irowp, &rows);
    3223:     ISInvertPermutation(colp, 0, &icolp);
    3224:     ISGetIndices(icolp, &cols);
    3225:     PetscMalloc(N * sizeof(PetscInt),         &cnew);
    3226:     PetscMalloc(N * sizeof(PetscScalar), &vnew);

    3228:     /* Setup bandwidth to include */
    3229:     if (band == PETSC_DECIDE) {
    3230:       if (frac <= 0.0)
    3231:         bw = (PetscInt) (M * 0.05);
    3232:       else
    3233:         bw = (PetscInt) (M * frac);
    3234:     } else {
    3235:       if (band <= 0) SETERRQ(PETSC_ERR_ARG_WRONG, "Bandwidth must be a positive integer");
    3236:       bw = band;
    3237:     }

    3239:     /* Put values into new matrix */
    3240:     MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, B);
    3241:     for(row = locRowStart, locRow = 0; row < locRowEnd; row++, locRow++) {
    3242:       MatGetRow(A, row, &nz, &cwork, &vwork);
    3243:       newRow   = rows[locRow]+locRowStart;
    3244:       for(col = 0, newNz = 0; col < nz; col++) {
    3245:         newCol = cols[cwork[col]];
    3246:         if ((newCol >= newRow - bw) && (newCol < newRow + bw) && (PetscAbsScalar(vwork[col]) >= tol)) {
    3247:           cnew[newNz] = newCol;
    3248:           vnew[newNz] = vwork[col];
    3249:           newNz++;
    3250:         }
    3251:       }
    3252:       MatSetValues(*B, 1, &newRow, newNz, cnew, vnew, INSERT_VALUES);
    3253:       MatRestoreRow(A, row, &nz, &cwork, &vwork);
    3254:     }
    3255:     PetscFree(cnew);
    3256:     PetscFree(vnew);
    3257:     MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
    3258:     MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
    3259:     ISRestoreIndices(irowp, &rows);
    3260:     ISRestoreIndices(icolp, &cols);
    3261:     ISDestroy(irowp);
    3262:     ISDestroy(icolp);
    3263:   } else {
    3264:     (*A->ops->permutesparsify)(A, band, frac, tol, rowp, colp, B);
    3265:   }
    3266:   PetscObjectStateIncrease((PetscObject)*B);
    3267:   return(0);
    3268: }

    3272: /*@
    3273:    MatEqual - Compares two matrices.

    3275:    Collective on Mat

    3277:    Input Parameters:
    3278: +  A - the first matrix
    3279: -  B - the second matrix

    3281:    Output Parameter:
    3282: .  flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise.

    3284:    Level: intermediate

    3286:    Concepts: matrices^equality between
    3287: @*/
    3288: PetscErrorCode PETSCMAT_DLLEXPORT MatEqual(Mat A,Mat B,PetscTruth *flg)
    3289: {

    3297:   MatPreallocated(B);
    3300:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    3301:   if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    3302:   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);
    3303:   if (!A->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",A->type_name);
    3304:   if (!B->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",B->type_name);
    3305:   if (A->ops->equal != B->ops->equal) SETERRQ2(PETSC_ERR_ARG_INCOMP,"A is type: %s\nB is type: %s",A->type_name,B->type_name);
    3306:   MatPreallocated(A);

    3308:   (*A->ops->equal)(A,B,flg);
    3309:   return(0);
    3310: }

    3314: /*@
    3315:    MatDiagonalScale - Scales a matrix on the left and right by diagonal
    3316:    matrices that are stored as vectors.  Either of the two scaling
    3317:    matrices can be PETSC_NULL.

    3319:    Collective on Mat

    3321:    Input Parameters:
    3322: +  mat - the matrix to be scaled
    3323: .  l - the left scaling vector (or PETSC_NULL)
    3324: -  r - the right scaling vector (or PETSC_NULL)

    3326:    Notes:
    3327:    MatDiagonalScale() computes A = LAR, where
    3328:    L = a diagonal matrix (stored as a vector), R = a diagonal matrix (stored as a vector)

    3330:    Level: intermediate

    3332:    Concepts: matrices^diagonal scaling
    3333:    Concepts: diagonal scaling of matrices

    3335: .seealso: MatScale()
    3336: @*/
    3337: PetscErrorCode PETSCMAT_DLLEXPORT MatDiagonalScale(Mat mat,Vec l,Vec r)
    3338: {

    3344:   if (!mat->ops->diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    3347:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    3348:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    3349:   MatPreallocated(mat);

    3351:   PetscLogEventBegin(MAT_Scale,mat,0,0,0);
    3352:   (*mat->ops->diagonalscale)(mat,l,r);
    3353:   PetscLogEventEnd(MAT_Scale,mat,0,0,0);
    3354:   PetscObjectStateIncrease((PetscObject)mat);
    3355:   return(0);
    3356: }

    3360: /*@
    3361:     MatScale - Scales all elements of a matrix by a given number.

    3363:     Collective on Mat

    3365:     Input Parameters:
    3366: +   mat - the matrix to be scaled
    3367: -   a  - the scaling value

    3369:     Output Parameter:
    3370: .   mat - the scaled matrix

    3372:     Level: intermediate

    3374:     Concepts: matrices^scaling all entries

    3376: .seealso: MatDiagonalScale()
    3377: @*/
    3378: PetscErrorCode PETSCMAT_DLLEXPORT MatScale(Mat mat,PetscScalar a)
    3379: {

    3385:   if (!mat->ops->scale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    3386:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    3387:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    3388:   MatPreallocated(mat);

    3390:   PetscLogEventBegin(MAT_Scale,mat,0,0,0);
    3391:   (*mat->ops->scale)(mat,a);
    3392:   PetscLogEventEnd(MAT_Scale,mat,0,0,0);
    3393:   PetscObjectStateIncrease((PetscObject)mat);
    3394:   return(0);
    3395: }

    3399: /*@
    3400:    MatNorm - Calculates various norms of a matrix.

    3402:    Collective on Mat

    3404:    Input Parameters:
    3405: +  mat - the matrix
    3406: -  type - the type of norm, NORM_1, NORM_FROBENIUS, NORM_INFINITY

    3408:    Output Parameters:
    3409: .  nrm - the resulting norm

    3411:    Level: intermediate

    3413:    Concepts: matrices^norm
    3414:    Concepts: norm^of matrix
    3415: @*/
    3416: PetscErrorCode PETSCMAT_DLLEXPORT MatNorm(Mat mat,NormType type,PetscReal *nrm)
    3417: {

    3425:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    3426:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    3427:   if (!mat->ops->norm) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    3428:   MatPreallocated(mat);

    3430:   (*mat->ops->norm)(mat,type,nrm);
    3431:   return(0);
    3432: }

    3434: /*
    3435:      This variable is used to prevent counting of MatAssemblyBegin() that
    3436:    are called from within a MatAssemblyEnd().
    3437: */
    3438: static PetscInt MatAssemblyEnd_InUse = 0;
    3441: /*@
    3442:    MatAssemblyBegin - Begins assembling the matrix.  This routine should
    3443:    be called after completing all calls to MatSetValues().

    3445:    Collective on Mat

    3447:    Input Parameters:
    3448: +  mat - the matrix
    3449: -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY
    3450:
    3451:    Notes:
    3452:    MatSetValues() generally caches the values.  The matrix is ready to
    3453:    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
    3454:    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
    3455:    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
    3456:    using the matrix.

    3458:    Level: beginner

    3460:    Concepts: matrices^assembling

    3462: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled()
    3463: @*/
    3464: PetscErrorCode PETSCMAT_DLLEXPORT MatAssemblyBegin(Mat mat,MatAssemblyType type)
    3465: {

    3471:   MatPreallocated(mat);
    3472:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.\nDid you forget to call MatSetUnfactored()?");
    3473:   if (mat->assembled) {
    3474:     mat->was_assembled = PETSC_TRUE;
    3475:     mat->assembled     = PETSC_FALSE;
    3476:   }
    3477:   if (!MatAssemblyEnd_InUse) {
    3478:     PetscLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);
    3479:     if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
    3480:     PetscLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);
    3481:   } else {
    3482:     if (mat->ops->assemblybegin){(*mat->ops->assemblybegin)(mat,type);}
    3483:   }
    3484:   return(0);
    3485: }

    3489: /*@
    3490:    MatAssembled - Indicates if a matrix has been assembled and is ready for
    3491:      use; for example, in matrix-vector product.

    3493:    Collective on Mat

    3495:    Input Parameter:
    3496: .  mat - the matrix

    3498:    Output Parameter:
    3499: .  assembled - PETSC_TRUE or PETSC_FALSE

    3501:    Level: advanced

    3503:    Concepts: matrices^assembled?

    3505: .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin()
    3506: @*/
    3507: PetscErrorCode PETSCMAT_DLLEXPORT MatAssembled(Mat mat,PetscTruth *assembled)
    3508: {
    3513:   *assembled = mat->assembled;
    3514:   return(0);
    3515: }

    3519: /*
    3520:     Processes command line options to determine if/how a matrix
    3521:   is to be viewed. Called by MatAssemblyEnd() and MatLoad().
    3522: */
    3523: PetscErrorCode MatView_Private(Mat mat)
    3524: {
    3525:   PetscErrorCode    ierr;
    3526:   PetscTruth        flg;
    3527:   static PetscTruth incall = PETSC_FALSE;

    3530:   if (incall) return(0);
    3531:   incall = PETSC_TRUE;
    3532:   PetscOptionsBegin(mat->comm,mat->prefix,"Matrix Options","Mat");
    3533:     PetscOptionsName("-mat_view_info","Information on matrix size","MatView",&flg);
    3534:     if (flg) {
    3535:       PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO);
    3536:       MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
    3537:       PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));
    3538:     }
    3539:     PetscOptionsName("-mat_view_info_detailed","Nonzeros in the matrix","MatView",&flg);
    3540:     if (flg) {
    3541:       PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO_DETAIL);
    3542:       MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
    3543:       PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));
    3544:     }
    3545:     PetscOptionsName("-mat_view","Print matrix to stdout","MatView",&flg);
    3546:     if (flg) {
    3547:       MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
    3548:     }
    3549:     PetscOptionsName("-mat_view_matlab","Print matrix to stdout in a format Matlab can read","MatView",&flg);
    3550:     if (flg) {
    3551:       PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_MATLAB);
    3552:       MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));
    3553:       PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));
    3554:     }
    3555: #if defined(PETSC_USE_SOCKET_VIEWER)
    3556:     PetscOptionsName("-mat_view_socket","Send matrix to socket (can be read from matlab)","MatView",&flg);
    3557:     if (flg) {
    3558:       MatView(mat,PETSC_VIEWER_SOCKET_(mat->comm));
    3559:       PetscViewerFlush(PETSC_VIEWER_SOCKET_(mat->comm));
    3560:     }
    3561: #endif
    3562:     PetscOptionsName("-mat_view_binary","Save matrix to file in binary format","MatView",&flg);
    3563:     if (flg) {
    3564:       MatView(mat,PETSC_VIEWER_BINARY_(mat->comm));
    3565:       PetscViewerFlush(PETSC_VIEWER_BINARY_(mat->comm));
    3566:     }
    3567:   PetscOptionsEnd();
    3568:   /* cannot have inside PetscOptionsBegin() because uses PetscOptionsBegin() */
    3569:   PetscOptionsHasName(mat->prefix,"-mat_view_draw",&flg);
    3570:   if (flg) {
    3571:     PetscOptionsHasName(mat->prefix,"-mat_view_contour",&flg);
    3572:     if (flg) {
    3573:       PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);
    3574:     }
    3575:     MatView(mat,PETSC_VIEWER_DRAW_(mat->comm));
    3576:     PetscViewerFlush(PETSC_VIEWER_DRAW_(mat->comm));
    3577:     if (flg) {
    3578:       PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));
    3579:     }
    3580:   }
    3581:   incall = PETSC_FALSE;
    3582:   return(0);
    3583: }

    3587: /*@
    3588:    MatAssemblyEnd - Completes assembling the matrix.  This routine should
    3589:    be called after MatAssemblyBegin().

    3591:    Collective on Mat

    3593:    Input Parameters:
    3594: +  mat - the matrix
    3595: -  type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY

    3597:    Options Database Keys:
    3598: +  -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly()
    3599: .  -mat_view_info_detailed - Prints more detailed info
    3600: .  -mat_view - Prints matrix in ASCII format
    3601: .  -mat_view_matlab - Prints matrix in Matlab format
    3602: .  -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX().
    3603: .  -display <name> - Sets display name (default is host)
    3604: .  -draw_pause <sec> - Sets number of seconds to pause after display
    3605: .  -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual)
    3606: .  -viewer_socket_machine <machine>
    3607: .  -viewer_socket_port <port>
    3608: .  -mat_view_binary - save matrix to file in binary format
    3609: -  -viewer_binary_filename <name>

    3611:    Notes:
    3612:    MatSetValues() generally caches the values.  The matrix is ready to
    3613:    use only after MatAssemblyBegin() and MatAssemblyEnd() have been called.
    3614:    Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES
    3615:    in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before
    3616:    using the matrix.

    3618:    Level: beginner

    3620: .seealso: MatAssemblyBegin(), MatSetValues(), PetscDrawOpenX(), MatView(), MatAssembled(), PetscViewerSocketOpen()
    3621: @*/
    3622: PetscErrorCode PETSCMAT_DLLEXPORT MatAssemblyEnd(Mat mat,MatAssemblyType type)
    3623: {
    3624:   PetscErrorCode  ierr;
    3625:   static PetscInt inassm = 0;
    3626:   PetscTruth      flg;

    3632:   inassm++;
    3633:   MatAssemblyEnd_InUse++;
    3634:   if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */
    3635:     PetscLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);
    3636:     if (mat->ops->assemblyend) {
    3637:       (*mat->ops->assemblyend)(mat,type);
    3638:     }
    3639:     PetscLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);
    3640:   } else {
    3641:     if (mat->ops->assemblyend) {
    3642:       (*mat->ops->assemblyend)(mat,type);
    3643:     }
    3644:   }

    3646:   /* Flush assembly is not a true assembly */
    3647:   if (type != MAT_FLUSH_ASSEMBLY) {
    3648:     mat->assembled  = PETSC_TRUE; mat->num_ass++;
    3649:   }
    3650:   mat->insertmode = NOT_SET_VALUES;
    3651:   MatAssemblyEnd_InUse--;
    3652:   PetscObjectStateIncrease((PetscObject)mat);
    3653:   if (!mat->symmetric_eternal) {
    3654:     mat->symmetric_set              = PETSC_FALSE;
    3655:     mat->hermitian_set              = PETSC_FALSE;
    3656:     mat->structurally_symmetric_set = PETSC_FALSE;
    3657:   }
    3658:   if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) {
    3659:     MatView_Private(mat);
    3660:     PetscOptionsHasName(mat->prefix,"-mat_is_symmetric",&flg);
    3661:     if (flg) {
    3662:       PetscReal tol = 0.0;
    3663:       PetscOptionsGetReal(mat->prefix,"-mat_is_symmetric",&tol,PETSC_NULL);
    3664:       MatIsSymmetric(mat,tol,&flg);
    3665:       if (flg) {
    3666:         PetscPrintf(mat->comm,"Matrix is symmetric (tolerance %g)\n",tol);
    3667:       } else {
    3668:         PetscPrintf(mat->comm,"Matrix is not symmetric (tolerance %g)\n",tol);
    3669:       }
    3670:     }
    3671:   }
    3672:   inassm--;
    3673:   PetscOptionsHasName(mat->prefix,"-help",&flg);
    3674:   if (flg) {
    3675:     MatPrintHelp(mat);
    3676:   }
    3677:   return(0);
    3678: }

    3683: /*@
    3684:    MatCompress - Tries to store the matrix in as little space as
    3685:    possible.  May fail if memory is already fully used, since it
    3686:    tries to allocate new space.

    3688:    Collective on Mat

    3690:    Input Parameters:
    3691: .  mat - the matrix

    3693:    Level: advanced

    3695: @*/
    3696: PetscErrorCode PETSCMAT_DLLEXPORT MatCompress(Mat mat)
    3697: {

    3703:   MatPreallocated(mat);
    3704:   if (mat->ops->compress) {(*mat->ops->compress)(mat);}
    3705:   return(0);
    3706: }

    3710: /*@
    3711:    MatSetOption - Sets a parameter option for a matrix. Some options
    3712:    may be specific to certain storage formats.  Some options
    3713:    determine how values will be inserted (or added). Sorted,
    3714:    row-oriented input will generally assemble the fastest. The default
    3715:    is row-oriented, nonsorted input.

    3717:    Collective on Mat

    3719:    Input Parameters:
    3720: +  mat - the matrix
    3721: -  option - the option, one of those listed below (and possibly others),
    3722:              e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR

    3724:    Options Describing Matrix Structure:
    3725: +    MAT_SYMMETRIC - symmetric in terms of both structure and value
    3726: .    MAT_HERMITIAN - transpose is the complex conjugation
    3727: .    MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure
    3728: .    MAT_NOT_SYMMETRIC - not symmetric in value
    3729: .    MAT_NOT_HERMITIAN - transpose is not the complex conjugation
    3730: .    MAT_NOT_STRUCTURALLY_SYMMETRIC - not symmetric nonzero structure
    3731: .    MAT_SYMMETRY_ETERNAL - if you would like the symmetry/Hermitian flag
    3732:                             you set to be kept with all future use of the matrix
    3733:                             including after MatAssemblyBegin/End() which could
    3734:                             potentially change the symmetry structure, i.e. you
    3735:                             KNOW the matrix will ALWAYS have the property you set.
    3736: -    MAT_NOT_SYMMETRY_ETERNAL - if MatAssemblyBegin/End() is called then the
    3737:                                 flags you set will be dropped (in case potentially
    3738:                                 the symmetry etc was lost).

    3740:    Options For Use with MatSetValues():
    3741:    Insert a logically dense subblock, which can be
    3742: +    MAT_ROW_ORIENTED - row-oriented (default)
    3743: .    MAT_COLUMN_ORIENTED - column-oriented
    3744: .    MAT_ROWS_SORTED - sorted by row
    3745: .    MAT_ROWS_UNSORTED - not sorted by row (default)
    3746: .    MAT_COLUMNS_SORTED - sorted by column
    3747: -    MAT_COLUMNS_UNSORTED - not sorted by column (default)

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

    3753:    When (re)assembling a matrix, we can restrict the input for
    3754:    efficiency/debugging purposes.  These options include
    3755: +    MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be
    3756:         allowed if they generate a new nonzero
    3757: .    MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed
    3758: .    MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if
    3759:          they generate a nonzero in a new diagonal (for block diagonal format only)
    3760: .    MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only)
    3761: .    MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries
    3762: .    MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry
    3763: -    MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly

    3765:    Notes:
    3766:    Some options are relevant only for particular matrix types and
    3767:    are thus ignored by others.  Other options are not supported by
    3768:    certain matrix types and will generate an error message if set.

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

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

    3781:    MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion
    3782:    that would generate a new entry in the nonzero structure instead produces
    3783:    an error. (Currently supported for AIJ and BAIJ formats only.)
    3784:    This is a useful flag when using SAME_NONZERO_PATTERN in calling
    3785:    KSPSetOperators() to ensure that the nonzero pattern truely does
    3786:    remain unchanged. Set after the first MatAssemblyEnd()

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

    3793:    MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for
    3794:    other processors should be dropped, rather than stashed.
    3795:    This is useful if you know that the "owning" processor is also
    3796:    always generating the correct matrix entries, so that PETSc need
    3797:    not transfer duplicate entries generated on another processor.
    3798:
    3799:    MAT_USE_HASH_TABLE indicates that a hash table be used to improve the
    3800:    searches during matrix assembly. When this flag is set, the hash table
    3801:    is created during the first Matrix Assembly. This hash table is
    3802:    used the next time through, during MatSetVaules()/MatSetVaulesBlocked()
    3803:    to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag
    3804:    should be used with MAT_USE_HASH_TABLE flag. This option is currently
    3805:    supported by MATMPIBAIJ format only.

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

    3810:    MAT_IGNORE_ZERO_ENTRIES - for AIJ/IS matrices this will stop zero values from creating
    3811:    a zero location in the matrix

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

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

    3819:    Level: intermediate

    3821:    Concepts: matrices^setting options

    3823: @*/
    3824: PetscErrorCode PETSCMAT_DLLEXPORT MatSetOption(Mat mat,MatOption op)
    3825: {

    3831:   MatPreallocated(mat);
    3832:   switch (op) {
    3833:   case MAT_SYMMETRIC:
    3834:     mat->symmetric                  = PETSC_TRUE;
    3835:     mat->structurally_symmetric     = PETSC_TRUE;
    3836:     mat->symmetric_set              = PETSC_TRUE;
    3837:     mat->structurally_symmetric_set = PETSC_TRUE;
    3838:     break;
    3839:   case MAT_HERMITIAN:
    3840:     mat->hermitian                  = PETSC_TRUE;
    3841:     mat->structurally_symmetric     = PETSC_TRUE;
    3842:     mat->hermitian_set              = PETSC_TRUE;
    3843:     mat->structurally_symmetric_set = PETSC_TRUE;
    3844:     break;
    3845:   case MAT_STRUCTURALLY_SYMMETRIC:
    3846:     mat->structurally_symmetric     = PETSC_TRUE;
    3847:     mat->structurally_symmetric_set = PETSC_TRUE;
    3848:     break;
    3849:   case MAT_NOT_SYMMETRIC:
    3850:     mat->symmetric                  = PETSC_FALSE;
    3851:     mat->symmetric_set              = PETSC_TRUE;
    3852:     break;
    3853:   case MAT_NOT_HERMITIAN:
    3854:     mat->hermitian                  = PETSC_FALSE;
    3855:     mat->hermitian_set              = PETSC_TRUE;
    3856:     break;
    3857:   case MAT_NOT_STRUCTURALLY_SYMMETRIC:
    3858:     mat->structurally_symmetric     = PETSC_FALSE;
    3859:     mat->structurally_symmetric_set = PETSC_TRUE;
    3860:     break;
    3861:   case MAT_SYMMETRY_ETERNAL:
    3862:     mat->symmetric_eternal          = PETSC_TRUE;
    3863:     break;
    3864:   case MAT_NOT_SYMMETRY_ETERNAL:
    3865:     mat->symmetric_eternal          = PETSC_FALSE;
    3866:     break;
    3867:   default:
    3868:     break;
    3869:   }
    3870:   if (mat->ops->setoption) {
    3871:     (*mat->ops->setoption)(mat,op);
    3872:   }
    3873:   return(0);
    3874: }

    3878: /*@
    3879:    MatZeroEntries - Zeros all entries of a matrix.  For sparse matrices
    3880:    this routine retains the old nonzero structure.

    3882:    Collective on Mat

    3884:    Input Parameters:
    3885: .  mat - the matrix

    3887:    Level: intermediate

    3889:    Concepts: matrices^zeroing

    3891: .seealso: MatZeroRows()
    3892: @*/
    3893: PetscErrorCode PETSCMAT_DLLEXPORT MatZeroEntries(Mat mat)
    3894: {

    3900:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    3901:   if (mat->insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for matrices where you have set values but not yet assembled");
    3902:   if (!mat->ops->zeroentries) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    3903:   MatPreallocated(mat);

    3905:   PetscLogEventBegin(MAT_ZeroEntries,mat,0,0,0);
    3906:   (*mat->ops->zeroentries)(mat);
    3907:   PetscLogEventEnd(MAT_ZeroEntries,mat,0,0,0);
    3908:   PetscObjectStateIncrease((PetscObject)mat);
    3909:   return(0);
    3910: }

    3914: /*@C
    3915:    MatZeroRows - Zeros all entries (except possibly the main diagonal)
    3916:    of a set of rows of a matrix.

    3918:    Collective on Mat

    3920:    Input Parameters:
    3921: +  mat - the matrix
    3922: .  numRows - the number of rows to remove
    3923: .  rows - the global row indices
    3924: -  diag - value put in all diagonals of eliminated rows

    3926:    Notes:
    3927:    For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
    3928:    but does not release memory.  For the dense and block diagonal
    3929:    formats this does not alter the nonzero structure.

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

    3935:    The user can set a value in the diagonal entry (or for the AIJ and
    3936:    row formats can optionally remove the main diagonal entry from the
    3937:    nonzero structure as well, by passing 0.0 as the final argument).

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

    3944:    Each processor should list the rows that IT wants zeroed

    3946:    Level: intermediate

    3948:    Concepts: matrices^zeroing rows

    3950: .seealso: MatZeroRowsIS(), MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
    3951: @*/
    3952: PetscErrorCode PETSCMAT_DLLEXPORT MatZeroRows(Mat mat,PetscInt numRows,const PetscInt rows[],PetscScalar diag)
    3953: {

    3960:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    3961:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    3962:   if (!mat->ops->zerorows) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    3963:   MatPreallocated(mat);

    3965:   (*mat->ops->zerorows)(mat,numRows,rows,diag);
    3966:   MatView_Private(mat);
    3967:   PetscObjectStateIncrease((PetscObject)mat);
    3968:   return(0);
    3969: }

    3973: /*@C
    3974:    MatZeroRowsIS - Zeros all entries (except possibly the main diagonal)
    3975:    of a set of rows of a matrix.

    3977:    Collective on Mat

    3979:    Input Parameters:
    3980: +  mat - the matrix
    3981: .  is - index set of rows to remove
    3982: -  diag - value put in all diagonals of eliminated rows

    3984:    Notes:
    3985:    For the AIJ and BAIJ matrix formats this removes the old nonzero structure,
    3986:    but does not release memory.  For the dense and block diagonal
    3987:    formats this does not alter the nonzero structure.

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

    3993:    The user can set a value in the diagonal entry (or for the AIJ and
    3994:    row formats can optionally remove the main diagonal entry from the
    3995:    nonzero structure as well, by passing 0.0 as the final argument).

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

    4002:    Each processor should list the rows that IT wants zeroed

    4004:    Level: intermediate

    4006:    Concepts: matrices^zeroing rows

    4008: .seealso: MatZeroRows(), MatZeroEntries(), MatZeroRowsLocal(), MatSetOption()
    4009: @*/
    4010: PetscErrorCode PETSCMAT_DLLEXPORT MatZeroRowsIS(Mat mat,IS is,PetscScalar diag)
    4011: {
    4012:   PetscInt       numRows;
    4013:   PetscInt       *rows;

    4020:   ISGetLocalSize(is,&numRows);
    4021:   ISGetIndices(is,&rows);
    4022:   MatZeroRows(mat,numRows,rows,diag);
    4023:   ISRestoreIndices(is,&rows);
    4024:   return(0);
    4025: }

    4029: /*@C
    4030:    MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
    4031:    of a set of rows of a matrix; using local numbering of rows.

    4033:    Collective on Mat

    4035:    Input Parameters:
    4036: +  mat - the matrix
    4037: .  numRows - the number of rows to remove
    4038: .  rows - the global row indices
    4039: -  diag - value put in all diagonals of eliminated rows

    4041:    Notes:
    4042:    Before calling MatZeroRowsLocal(), the user must first set the
    4043:    local-to-global mapping by calling MatSetLocalToGlobalMapping().

    4045:    For the AIJ matrix formats this removes the old nonzero structure,
    4046:    but does not release memory.  For the dense and block diagonal
    4047:    formats this does not alter the nonzero structure.

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

    4053:    The user can set a value in the diagonal entry (or for the AIJ and
    4054:    row formats can optionally remove the main diagonal entry from the
    4055:    nonzero structure as well, by passing 0.0 as the final argument).

    4057:    Level: intermediate

    4059:    Concepts: matrices^zeroing

    4061: .seealso: MatZeroRows(), MatZeroRowsLocalIS(), MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
    4062: @*/
    4063: PetscErrorCode PETSCMAT_DLLEXPORT MatZeroRowsLocal(Mat mat,PetscInt numRows,const PetscInt rows[],PetscScalar diag)
    4064: {

    4071:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    4072:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    4073:   MatPreallocated(mat);

    4075:   if (mat->ops->zerorowslocal) {
    4076:     (*mat->ops->zerorowslocal)(mat,numRows,rows,diag);
    4077:   } else {
    4078:     IS is, newis;
    4079:     PetscInt *newRows;

    4081:     if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Need to provide local to global mapping to matrix first");
    4082:     ISCreateGeneral(PETSC_COMM_SELF,numRows,rows,&is);
    4083:     ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);
    4084:     ISGetIndices(newis,&newRows);
    4085:     (*mat->ops->zerorows)(mat,numRows,newRows,diag);
    4086:     ISRestoreIndices(newis,&newRows);
    4087:     ISDestroy(newis);
    4088:     ISDestroy(is);
    4089:   }
    4090:   PetscObjectStateIncrease((PetscObject)mat);
    4091:   return(0);
    4092: }

    4096: /*@C
    4097:    MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal)
    4098:    of a set of rows of a matrix; using local numbering of rows.

    4100:    Collective on Mat

    4102:    Input Parameters:
    4103: +  mat - the matrix
    4104: .  is - index set of rows to remove
    4105: -  diag - value put in all diagonals of eliminated rows

    4107:    Notes:
    4108:    Before calling MatZeroRowsLocal(), the user must first set the
    4109:    local-to-global mapping by calling MatSetLocalToGlobalMapping().

    4111:    For the AIJ matrix formats this removes the old nonzero structure,
    4112:    but does not release memory.  For the dense and block diagonal
    4113:    formats this does not alter the nonzero structure.

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

    4119:    The user can set a value in the diagonal entry (or for the AIJ and
    4120:    row formats can optionally remove the main diagonal entry from the
    4121:    nonzero structure as well, by passing 0.0 as the final argument).

    4123:    Level: intermediate

    4125:    Concepts: matrices^zeroing

    4127: .seealso: MatZeroRows(), MatZeroRowsLocal(), MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping
    4128: @*/
    4129: PetscErrorCode PETSCMAT_DLLEXPORT MatZeroRowsLocalIS(Mat mat,IS is,PetscScalar diag)
    4130: {
    4132:   PetscInt       numRows;
    4133:   PetscInt       *rows;

    4139:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    4140:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    4141:   MatPreallocated(mat);

    4143:   ISGetLocalSize(is,&numRows);
    4144:   ISGetIndices(is,&rows);
    4145:   MatZeroRowsLocal(mat,numRows,rows,diag);
    4146:   ISRestoreIndices(is,&rows);
    4147:   return(0);
    4148: }

    4152: /*@
    4153:    MatGetSize - Returns the numbers of rows and columns in a matrix.

    4155:    Not Collective

    4157:    Input Parameter:
    4158: .  mat - the matrix

    4160:    Output Parameters:
    4161: +  m - the number of global rows
    4162: -  n - the number of global columns

    4164:    Note: both output parameters can be PETSC_NULL on input.

    4166:    Level: beginner

    4168:    Concepts: matrices^size

    4170: .seealso: MatGetLocalSize()
    4171: @*/
    4172: PetscErrorCode PETSCMAT_DLLEXPORT MatGetSize(Mat mat,PetscInt *m,PetscInt* n)
    4173: {
    4176:   if (m) *m = mat->M;
    4177:   if (n) *n = mat->N;
    4178:   return(0);
    4179: }

    4183: /*@
    4184:    MatGetLocalSize - Returns the number of rows and columns in a matrix
    4185:    stored locally.  This information may be implementation dependent, so
    4186:    use with care.

    4188:    Not Collective

    4190:    Input Parameters:
    4191: .  mat - the matrix

    4193:    Output Parameters:
    4194: +  m - the number of local rows
    4195: -  n - the number of local columns

    4197:    Note: both output parameters can be PETSC_NULL on input.

    4199:    Level: beginner

    4201:    Concepts: matrices^local size

    4203: .seealso: MatGetSize()
    4204: @*/
    4205: PetscErrorCode PETSCMAT_DLLEXPORT MatGetLocalSize(Mat mat,PetscInt *m,PetscInt* n)
    4206: {
    4211:   if (m) *m = mat->m;
    4212:   if (n) *n = mat->n;
    4213:   return(0);
    4214: }

    4218: /*@
    4219:    MatGetOwnershipRange - Returns the range of matrix rows owned by
    4220:    this processor, assuming that the matrix is laid out with the first
    4221:    n1 rows on the first processor, the next n2 rows on the second, etc.
    4222:    For certain parallel layouts this range may not be well defined.

    4224:    Not Collective

    4226:    Input Parameters:
    4227: .  mat - the matrix

    4229:    Output Parameters:
    4230: +  m - the global index of the first local row
    4231: -  n - one more than the global index of the last local row

    4233:    Note: both output parameters can be PETSC_NULL on input.

    4235:    Level: beginner

    4237:    Concepts: matrices^row ownership
    4238: @*/
    4239: PetscErrorCode PETSCMAT_DLLEXPORT MatGetOwnershipRange(Mat mat,PetscInt *m,PetscInt* n)
    4240: {

    4248:   MatPreallocated(mat);
    4249:   PetscMapGetLocalRange(mat->rmap,m,n);
    4250:   return(0);
    4251: }

    4255: /*@
    4256:    MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix.
    4257:    Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric()
    4258:    to complete the factorization.

    4260:    Collective on Mat

    4262:    Input Parameters:
    4263: +  mat - the matrix
    4264: .  row - row permutation
    4265: .  column - column permutation
    4266: -  info - structure containing
    4267: $      levels - number of levels of fill.
    4268: $      expected fill - as ratio of original fill.
    4269: $      1 or 0 - indicating force fill on diagonal (improves robustness for matrices
    4270:                 missing diagonal entries)

    4272:    Output Parameters:
    4273: .  fact - new matrix that has been symbolically factored

    4275:    Notes:
    4276:    See the users manual for additional information about
    4277:    choosing the fill factor for better efficiency.

    4279:    Most users should employ the simplified KSP interface for linear solvers
    4280:    instead of working directly with matrix algebra routines such as this.
    4281:    See, e.g., KSPCreate().

    4283:    Level: developer

    4285:   Concepts: matrices^symbolic LU factorization
    4286:   Concepts: matrices^factorization
    4287:   Concepts: LU^symbolic factorization

    4289: .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
    4290:           MatGetOrdering(), MatFactorInfo

    4292: @*/
    4293: PetscErrorCode PETSCMAT_DLLEXPORT MatILUFactorSymbolic(Mat mat,IS row,IS col,MatFactorInfo *info,Mat *fact)
    4294: {

    4304:   if (info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels of fill negative %D",(PetscInt)info->levels);
    4305:   if (info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",info->fill);
    4306:   if (!mat->ops->ilufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s  symbolic ILU",mat->type_name);
    4307:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    4308:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    4309:   MatPreallocated(mat);

    4311:   PetscLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);
    4312:   (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);
    4313:   PetscLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);
    4314:   return(0);
    4315: }

    4319: /*@
    4320:    MatICCFactorSymbolic - Performs symbolic incomplete
    4321:    Cholesky factorization for a symmetric matrix.  Use
    4322:    MatCholeskyFactorNumeric() to complete the factorization.

    4324:    Collective on Mat

    4326:    Input Parameters:
    4327: +  mat - the matrix
    4328: .  perm - row and column permutation
    4329: -  info - structure containing
    4330: $      levels - number of levels of fill.
    4331: $      expected fill - as ratio of original fill.

    4333:    Output Parameter:
    4334: .  fact - the factored matrix

    4336:    Notes:
    4337:    Currently only no-fill factorization is supported.

    4339:    Most users should employ the simplified KSP interface for linear solvers
    4340:    instead of working directly with matrix algebra routines such as this.
    4341:    See, e.g., KSPCreate().

    4343:    Level: developer

    4345:   Concepts: matrices^symbolic incomplete Cholesky factorization
    4346:   Concepts: matrices^factorization
    4347:   Concepts: Cholsky^symbolic factorization

    4349: .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor(), MatFactorInfo
    4350: @*/
    4351: PetscErrorCode PETSCMAT_DLLEXPORT MatICCFactorSymbolic(Mat mat,IS perm,MatFactorInfo *info,Mat *fact)
    4352: {

    4361:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    4362:   if (info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels negative %D",(PetscInt) info->levels);
    4363:   if (info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",info->fill);
    4364:   if (!mat->ops->iccfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s  symbolic ICC",mat->type_name);
    4365:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    4366:   MatPreallocated(mat);

    4368:   PetscLogEventBegin(MAT_ICCFactorSymbolic,mat,perm,0,0);
    4369:   (*mat->ops->iccfactorsymbolic)(mat,perm,info,fact);
    4370:   PetscLogEventEnd(MAT_ICCFactorSymbolic,mat,perm,0,0);
    4371:   return(0);
    4372: }

    4376: /*@C
    4377:    MatGetArray - Returns a pointer to the element values in the matrix.
    4378:    The result of this routine is dependent on the underlying matrix data
    4379:    structure, and may not even work for certain matrix types.  You MUST
    4380:    call MatRestoreArray() when you no longer need to access the array.

    4382:    Not Collective

    4384:    Input Parameter:
    4385: .  mat - the matrix

    4387:    Output Parameter:
    4388: .  v - the location of the values

    4391:    Fortran Note:
    4392:    This routine is used differently from Fortran, e.g.,
    4393: .vb
    4394:         Mat         mat
    4395:         PetscScalar mat_array(1)
    4396:         PetscOffset i_mat
    4397:         PetscErrorCode ierr
    4398:         call MatGetArray(mat,mat_array,i_mat,ierr)

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

    4404:         [... other code ...]
    4405:         call MatRestoreArray(mat,mat_array,i_mat,ierr)
    4406: .ve

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

    4411:    Level: advanced

    4413:    Concepts: matrices^access array

    4415: .seealso: MatRestoreArray(), MatGetArrayF90()
    4416: @*/
    4417: PetscErrorCode PETSCMAT_DLLEXPORT MatGetArray(Mat mat,PetscScalar *v[])
    4418: {

    4425:   if (!mat->ops->getarray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    4426:   MatPreallocated(mat);
    4427:   (*mat->ops->getarray)(mat,v);
    4428:   return(0);
    4429: }

    4433: /*@C
    4434:    MatRestoreArray - Restores the matrix after MatGetArray() has been called.

    4436:    Not Collective

    4438:    Input Parameter:
    4439: +  mat - the matrix
    4440: -  v - the location of the values

    4442:    Fortran Note:
    4443:    This routine is used differently from Fortran, e.g.,
    4444: .vb
    4445:         Mat         mat
    4446:         PetscScalar mat_array(1)
    4447:         PetscOffset i_mat
    4448:         PetscErrorCode ierr
    4449:         call MatGetArray(mat,mat_array,i_mat,ierr)

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

    4455:         [... other code ...]
    4456:         call MatRestoreArray(mat,mat_array,i_mat,ierr)
    4457: .ve

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

    4462:    Level: advanced

    4464: .seealso: MatGetArray(), MatRestoreArrayF90()
    4465: @*/
    4466: PetscErrorCode PETSCMAT_DLLEXPORT MatRestoreArray(Mat mat,PetscScalar *v[])
    4467: {

    4474: #if defined(PETSC_USE_DEBUG)
    4475:   CHKMEMQ;
    4476: #endif
    4477:   if (!mat->ops->restorearray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    4478:   (*mat->ops->restorearray)(mat,v);
    4479:   PetscObjectStateIncrease((PetscObject)mat);
    4480:   return(0);
    4481: }

    4485: /*@C
    4486:    MatGetSubMatrices - Extracts several submatrices from a matrix. If submat
    4487:    points to an array of valid matrices, they may be reused to store the new
    4488:    submatrices.

    4490:    Collective on Mat

    4492:    Input Parameters:
    4493: +  mat - the matrix
    4494: .  n   - the number of submatrixes to be extracted (on this processor, may be zero)
    4495: .  irow, icol - index sets of rows and columns to extract
    4496: -  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

    4498:    Output Parameter:
    4499: .  submat - the array of submatrices

    4501:    Notes:
    4502:    MatGetSubMatrices() can extract only sequential submatrices
    4503:    (from both sequential and parallel matrices). Use MatGetSubMatrix()
    4504:    to extract a parallel submatrix.

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

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

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

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

    4519:    For BAIJ matrices the index sets must respect the block structure, that is if they
    4520:    request one row/column in a block, they must request all rows/columns that are in
    4521:    that block. For example, if the block size is 2 you cannot request just row 0 and
    4522:    column 0.

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

    4528:    Level: advanced

    4530:    Concepts: matrices^accessing submatrices
    4531:    Concepts: submatrices

    4533: .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal()
    4534: @*/
    4535: PetscErrorCode PETSCMAT_DLLEXPORT MatGetSubMatrices(Mat mat,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *submat[])
    4536: {
    4538:   PetscInt        i;
    4539:   PetscTruth      eq;

    4544:   if (n) {
    4549:   }
    4551:   if (n && scall == MAT_REUSE_MATRIX) {
    4554:   }
    4555:   if (!mat->ops->getsubmatrices) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    4556:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    4557:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    4558:   MatPreallocated(mat);

    4560:   PetscLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);
    4561:   (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);
    4562:   PetscLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);
    4563:   for (i=0; i<n; i++) {
    4564:     if (mat->symmetric || mat->structurally_symmetric || mat->hermitian) {
    4565:       ISEqual(irow[i],icol[i],&eq);
    4566:       if (eq) {
    4567:         if (mat->symmetric){
    4568:           MatSetOption((*submat)[i],MAT_SYMMETRIC);
    4569:         } else if (mat->hermitian) {
    4570:           MatSetOption((*submat)[i],MAT_HERMITIAN);
    4571:         } else if (mat->structurally_symmetric) {
    4572:           MatSetOption((*submat)[i],MAT_STRUCTURALLY_SYMMETRIC);
    4573:         }
    4574:       }
    4575:     }
    4576:   }
    4577:   return(0);
    4578: }

    4582: /*@C
    4583:    MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices().

    4585:    Collective on Mat

    4587:    Input Parameters:
    4588: +  n - the number of local matrices
    4589: -  mat - the matrices (note that this is a pointer to the array of matrices, just to match the calling
    4590:                        sequence of MatGetSubMatrices())

    4592:    Level: advanced

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

    4596: .seealso: MatGetSubMatrices()
    4597: @*/
    4598: PetscErrorCode PETSCMAT_DLLEXPORT MatDestroyMatrices(PetscInt n,Mat *mat[])
    4599: {
    4601:   PetscInt       i;

    4604:   if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of matrices %D",n);
    4606:   for (i=0; i<n; i++) {
    4607:     MatDestroy((*mat)[i]);
    4608:   }
    4609:   /* memory is allocated even if n = 0 */
    4610:   PetscFree(*mat);
    4611:   return(0);
    4612: }

    4616: /*@
    4617:    MatIncreaseOverlap - Given a set of submatrices indicated by index sets,
    4618:    replaces the index sets by larger ones that represent submatrices with
    4619:    additional overlap.

    4621:    Collective on Mat

    4623:    Input Parameters:
    4624: +  mat - the matrix
    4625: .  n   - the number of index sets
    4626: .  is  - the array of index sets (these index sets will changed during the call)
    4627: -  ov  - the additional overlap requested

    4629:    Level: developer

    4631:    Concepts: overlap
    4632:    Concepts: ASM^computing overlap

    4634: .seealso: MatGetSubMatrices()
    4635: @*/
    4636: PetscErrorCode PETSCMAT_DLLEXPORT MatIncreaseOverlap(Mat mat,PetscInt n,IS is[],PetscInt ov)
    4637: {

    4643:   if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more domains, you have %D",n);
    4644:   if (n) {
    4647:   }
    4648:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    4649:   if (mat->factor)     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    4650:   MatPreallocated(mat);

    4652:   if (!ov) return(0);
    4653:   if (!mat->ops->increaseoverlap) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    4654:   PetscLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);
    4655:   (*mat->ops->increaseoverlap)(mat,n,is,ov);
    4656:   PetscLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);
    4657:   return(0);
    4658: }

    4662: /*@
    4663:    MatPrintHelp - Prints all the options for the matrix.

    4665:    Collective on Mat

    4667:    Input Parameter:
    4668: .  mat - the matrix

    4670:    Options Database Keys:
    4671: +  -help - Prints matrix options
    4672: -  -h - Prints matrix options

    4674:    Level: developer

    4676: .seealso: MatCreate(), MatCreateXXX()
    4677: @*/
    4678: PetscErrorCode PETSCMAT_DLLEXPORT MatPrintHelp(Mat mat)
    4679: {
    4680:   static PetscTruth called = PETSC_FALSE;
    4681:   PetscErrorCode    ierr;

    4686:   MatPreallocated(mat);

    4688:   if (!called) {
    4689:     if (mat->ops->printhelp) {
    4690:       (*mat->ops->printhelp)(mat);
    4691:     }
    4692:     called = PETSC_TRUE;
    4693:   }
    4694:   return(0);
    4695: }

    4699: /*@
    4700:    MatGetBlockSize - Returns the matrix block size; useful especially for the
    4701:    block row and block diagonal formats.
    4702:
    4703:    Not Collective

    4705:    Input Parameter:
    4706: .  mat - the matrix

    4708:    Output Parameter:
    4709: .  bs - block size

    4711:    Notes:
    4712:    Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG.
    4713:    Block row formats are MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ

    4715:    Level: intermediate

    4717:    Concepts: matrices^block size

    4719: .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag()
    4720: @*/
    4721: PetscErrorCode PETSCMAT_DLLEXPORT MatGetBlockSize(Mat mat,PetscInt *bs)
    4722: {

    4729:   MatPreallocated(mat);
    4730:   *bs = mat->bs;
    4731:   return(0);
    4732: }

    4736: /*@
    4737:    MatSetBlockSize - Sets the matrix block size; for many matrix types you
    4738:      cannot use this and MUST set the blocksize when you preallocate the matrix
    4739:
    4740:    Not Collective

    4742:    Input Parameters:
    4743: +  mat - the matrix
    4744: -  bs - block size

    4746:    Notes:
    4747:      Only works for shell and AIJ matrices

    4749:    Level: intermediate

    4751:    Concepts: matrices^block size

    4753: .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag(), MatGetBlockSize()
    4754: @*/
    4755: PetscErrorCode PETSCMAT_DLLEXPORT MatSetBlockSize(Mat mat,PetscInt bs)
    4756: {

    4762:   MatPreallocated(mat);
    4763:   if (mat->ops->setblocksize) {
    4764:     mat->bs = bs;
    4765:     (*mat->ops->setblocksize)(mat,bs);
    4766:   } else {
    4767:     SETERRQ1(PETSC_ERR_ARG_INCOMP,"Cannot set the blocksize for matrix type %s",mat->type_name);
    4768:   }
    4769:   return(0);
    4770: }

    4774: /*@C
    4775:     MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices.

    4777:    Collective on Mat

    4779:     Input Parameters:
    4780: +   mat - the matrix
    4781: .   shift -  0 or 1 indicating we want the indices starting at 0 or 1
    4782: -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
    4783:                 symmetrized

    4785:     Output Parameters:
    4786: +   n - number of rows in the (possibly compressed) matrix
    4787: .   ia - the row pointers
    4788: .   ja - the column indices
    4789: -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned

    4791:     Level: developer

    4793: .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
    4794: @*/
    4795: PetscErrorCode PETSCMAT_DLLEXPORT MatGetRowIJ(Mat mat,PetscInt shift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt* ja[],PetscTruth *done)
    4796: {

    4806:   MatPreallocated(mat);
    4807:   if (!mat->ops->getrowij) *done = PETSC_FALSE;
    4808:   else {
    4809:     *done = PETSC_TRUE;
    4810:     (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);
    4811:   }
    4812:   return(0);
    4813: }

    4817: /*@C
    4818:     MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices.

    4820:     Collective on Mat

    4822:     Input Parameters:
    4823: +   mat - the matrix
    4824: .   shift - 1 or zero indicating we want the indices starting at 0 or 1
    4825: -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
    4826:                 symmetrized

    4828:     Output Parameters:
    4829: +   n - number of columns in the (possibly compressed) matrix
    4830: .   ia - the column pointers
    4831: .   ja - the row indices
    4832: -   done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned

    4834:     Level: developer

    4836: .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
    4837: @*/
    4838: PetscErrorCode PETSCMAT_DLLEXPORT MatGetColumnIJ(Mat mat,PetscInt shift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt* ja[],PetscTruth *done)
    4839: {

    4849:   MatPreallocated(mat);
    4850:   if (!mat->ops->getcolumnij) *done = PETSC_FALSE;
    4851:   else {
    4852:     *done = PETSC_TRUE;
    4853:     (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);
    4854:   }
    4855:   return(0);
    4856: }

    4860: /*@C
    4861:     MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with
    4862:     MatGetRowIJ().

    4864:     Collective on Mat

    4866:     Input Parameters:
    4867: +   mat - the matrix
    4868: .   shift - 1 or zero indicating we want the indices starting at 0 or 1
    4869: -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
    4870:                 symmetrized

    4872:     Output Parameters:
    4873: +   n - size of (possibly compressed) matrix
    4874: .   ia - the row pointers
    4875: .   ja - the column indices
    4876: -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned

    4878:     Level: developer

    4880: .seealso: MatGetRowIJ(), MatRestoreColumnIJ()
    4881: @*/
    4882: PetscErrorCode PETSCMAT_DLLEXPORT MatRestoreRowIJ(Mat mat,PetscInt shift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt* ja[],PetscTruth *done)
    4883: {

    4892:   MatPreallocated(mat);

    4894:   if (!mat->ops->restorerowij) *done = PETSC_FALSE;
    4895:   else {
    4896:     *done = PETSC_TRUE;
    4897:     (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);
    4898:   }
    4899:   return(0);
    4900: }

    4904: /*@C
    4905:     MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with
    4906:     MatGetColumnIJ().

    4908:     Collective on Mat

    4910:     Input Parameters:
    4911: +   mat - the matrix
    4912: .   shift - 1 or zero indicating we want the indices starting at 0 or 1
    4913: -   symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be
    4914:                 symmetrized

    4916:     Output Parameters:
    4917: +   n - size of (possibly compressed) matrix
    4918: .   ia - the column pointers
    4919: .   ja - the row indices
    4920: -   done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned

    4922:     Level: developer

    4924: .seealso: MatGetColumnIJ(), MatRestoreRowIJ()
    4925: @*/
    4926: PetscErrorCode PETSCMAT_DLLEXPORT MatRestoreColumnIJ(Mat mat,PetscInt shift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt* ja[],PetscTruth *done)
    4927: {

    4936:   MatPreallocated(mat);

    4938:   if (!mat->ops->restorecolumnij) *done = PETSC_FALSE;
    4939:   else {
    4940:     *done = PETSC_TRUE;
    4941:     (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);
    4942:   }
    4943:   return(0);
    4944: }

    4948: /*@C
    4949:     MatColoringPatch -Used inside matrix coloring routines that
    4950:     use MatGetRowIJ() and/or MatGetColumnIJ().

    4952:     Collective on Mat

    4954:     Input Parameters:
    4955: +   mat - the matrix
    4956: .   n   - number of colors
    4957: -   colorarray - array indicating color for each column

    4959:     Output Parameters:
    4960: .   iscoloring - coloring generated using colorarray information

    4962:     Level: developer

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

    4966: @*/
    4967: PetscErrorCode PETSCMAT_DLLEXPORT MatColoringPatch(Mat mat,PetscInt n,PetscInt ncolors,ISColoringValue colorarray[],ISColoring *iscoloring)
    4968: {

    4976:   MatPreallocated(mat);

    4978:   if (!mat->ops->coloringpatch){
    4979:     ISColoringCreate(mat->comm,n,colorarray,iscoloring);
    4980:   } else {
    4981:     (*mat->ops->coloringpatch)(mat,n,ncolors,colorarray,iscoloring);
    4982:   }
    4983:   return(0);
    4984: }

    4989: /*@
    4990:    MatSetUnfactored - Resets a factored matrix to be treated as unfactored.

    4992:    Collective on Mat

    4994:    Input Parameter:
    4995: .  mat - the factored matrix to be reset

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

    5004:    Note that one can specify in-place ILU(0) factorization by calling
    5005: .vb
    5006:      PCType(pc,PCILU);
    5007:      PCILUSeUseInPlace(pc);
    5008: .ve
    5009:    or by using the options -pc_type ilu -pc_ilu_in_place

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

    5017:    Most users should employ the simplified KSP interface for linear solvers
    5018:    instead of working directly with matrix algebra routines such as this.
    5019:    See, e.g., KSPCreate().

    5021:    Level: developer

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

    5025:    Concepts: matrices^unfactored

    5027: @*/
    5028: PetscErrorCode PETSCMAT_DLLEXPORT MatSetUnfactored(Mat mat)
    5029: {

    5035:   MatPreallocated(mat);
    5036:   mat->factor = 0;
    5037:   if (!mat->ops->setunfactored) return(0);
    5038:   (*mat->ops->setunfactored)(mat);
    5039:   return(0);
    5040: }

    5042: /*MC
    5043:     MatGetArrayF90 - Accesses a matrix array from Fortran90.

    5045:     Synopsis:
    5046:     MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)

    5048:     Not collective

    5050:     Input Parameter:
    5051: .   x - matrix

    5053:     Output Parameters:
    5054: +   xx_v - the Fortran90 pointer to the array
    5055: -   ierr - error code

    5057:     Example of Usage:
    5058: .vb
    5059:       PetscScalar, pointer xx_v(:)
    5060:       ....
    5061:       call MatGetArrayF90(x,xx_v,ierr)
    5062:       a = xx_v(3)
    5063:       call MatRestoreArrayF90(x,xx_v,ierr)
    5064: .ve

    5066:     Notes:
    5067:     Not yet supported for all F90 compilers

    5069:     Level: advanced

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

    5073:     Concepts: matrices^accessing array

    5075: M*/

    5077: /*MC
    5078:     MatRestoreArrayF90 - Restores a matrix array that has been
    5079:     accessed with MatGetArrayF90().

    5081:     Synopsis:
    5082:     MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr)

    5084:     Not collective

    5086:     Input Parameters:
    5087: +   x - matrix
    5088: -   xx_v - the Fortran90 pointer to the array

    5090:     Output Parameter:
    5091: .   ierr - error code

    5093:     Example of Usage:
    5094: .vb
    5095:        PetscScalar, pointer xx_v(:)
    5096:        ....
    5097:        call MatGetArrayF90(x,xx_v,ierr)
    5098:        a = xx_v(3)
    5099:        call MatRestoreArrayF90(x,xx_v,ierr)
    5100: .ve
    5101:
    5102:     Notes:
    5103:     Not yet supported for all F90 compilers

    5105:     Level: advanced

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

    5109: M*/

    5114: /*@
    5115:     MatGetSubMatrix - Gets a single submatrix on the same number of processors
    5116:                       as the original matrix.

    5118:     Collective on Mat

    5120:     Input Parameters:
    5121: +   mat - the original matrix
    5122: .   isrow - rows this processor should obtain
    5123: .   iscol - columns for all processors you wish to keep
    5124: .   csize - number of columns "local" to this processor (does nothing for sequential
    5125:             matrices). This should match the result from VecGetLocalSize(x,...) if you
    5126:             plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE
    5127: -   cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX

    5129:     Output Parameter:
    5130: .   newmat - the new submatrix, of the same type as the old

    5132:     Level: advanced

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

    5137:       The first time this is called you should use a cll of MAT_INITIAL_MATRIX,
    5138:    the MatGetSubMatrix() routine will create the newmat for you. Any additional calls
    5139:    to this routine with a mat of the same nonzero structure and with a cll of MAT_REUSE_MATRIX
    5140:    will reuse the matrix generated the first time.

    5142:     Concepts: matrices^submatrices

    5144: .seealso: MatGetSubMatrices(), ISAllGather()
    5145: @*/
    5146: PetscErrorCode PETSCMAT_DLLEXPORT MatGetSubMatrix(Mat mat,IS isrow,IS iscol,PetscInt csize,MatReuse cll,Mat *newmat)
    5147: {
    5149:   PetscMPIInt    size;
    5150:   Mat            *local;

    5159:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    5160:   MatPreallocated(mat);
    5161:   MPI_Comm_size(mat->comm,&size);

    5163:   /* if original matrix is on just one processor then use submatrix generated */
    5164:   if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) {
    5165:     MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);
    5166:     return(0);
    5167:   } else if (!mat->ops->getsubmatrix && size == 1) {
    5168:     MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);
    5169:     *newmat = *local;
    5170:     PetscFree(local);
    5171:     return(0);
    5172:   }

    5174:   if (!mat->ops->getsubmatrix) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    5175:   (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);
    5176:   PetscObjectStateIncrease((PetscObject)*newmat);
    5177:   return(0);
    5178: }

    5182: /*@C
    5183:    MatGetPetscMaps - Returns the maps associated with the matrix.

    5185:    Not Collective

    5187:    Input Parameter:
    5188: .  mat - the matrix

    5190:    Output Parameters:
    5191: +  rmap - the row (right) map
    5192: -  cmap - the column (left) map

    5194:    Level: developer

    5196:    Concepts: maps^getting from matrix

    5198: @*/
    5199: PetscErrorCode PETSCMAT_DLLEXPORT MatGetPetscMaps(Mat mat,PetscMap *rmap,PetscMap *cmap)
    5200: {

    5206:   MatPreallocated(mat);
    5207:   (*mat->ops->getmaps)(mat,rmap,cmap);
    5208:   return(0);
    5209: }

    5211: /*
    5212:       Version that works for all PETSc matrices
    5213: */
    5216: PetscErrorCode MatGetPetscMaps_Petsc(Mat mat,PetscMap *rmap,PetscMap *cmap)
    5217: {
    5219:   if (rmap) *rmap = mat->rmap;
    5220:   if (cmap) *cmap = mat->cmap;
    5221:   return(0);
    5222: }

    5226: /*@
    5227:    MatStashSetInitialSize - sets the sizes of the matrix stash, that is
    5228:    used during the assembly process to store values that belong to
    5229:    other processors.

    5231:    Not Collective

    5233:    Input Parameters:
    5234: +  mat   - the matrix
    5235: .  size  - the initial size of the stash.
    5236: -  bsize - the initial size of the block-stash(if used).

    5238:    Options Database Keys:
    5239: +   -matstash_initial_size <size> or <size0,size1,...sizep-1>
    5240: -   -matstash_block_initial_size <bsize>  or <bsize0,bsize1,...bsizep-1>

    5242:    Level: intermediate

    5244:    Notes:
    5245:      The block-stash is used for values set with MatSetValuesBlocked() while
    5246:      the stash is used for values set with MatSetValues()

    5248:      Run with the option -log_info and look for output of the form
    5249:      MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
    5250:      to determine the appropriate value, MM, to use for size and
    5251:      MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
    5252:      to determine the value, BMM to use for bsize

    5254:    Concepts: stash^setting matrix size
    5255:    Concepts: matrices^stash

    5257: @*/
    5258: PetscErrorCode PETSCMAT_DLLEXPORT MatStashSetInitialSize(Mat mat,PetscInt size, PetscInt bsize)
    5259: {

    5265:   MatStashSetInitialSize_Private(&mat->stash,size);
    5266:   MatStashSetInitialSize_Private(&mat->bstash,bsize);
    5267:   return(0);
    5268: }

    5272: /*@
    5273:    MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of
    5274:      the matrix

    5276:    Collective on Mat

    5278:    Input Parameters:
    5279: +  mat   - the matrix
    5280: .  x,y - the vectors
    5281: -  w - where the result is stored

    5283:    Level: intermediate

    5285:    Notes:
    5286:     w may be the same vector as y.

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

    5291:     Concepts: interpolation

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

    5295: @*/
    5296: PetscErrorCode PETSCMAT_DLLEXPORT MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w)
    5297: {
    5299:   PetscInt       M,N;

    5307:   MatPreallocated(A);
    5308:   MatGetSize(A,&M,&N);
    5309:   if (N > M) {
    5310:     MatMultTransposeAdd(A,x,y,w);
    5311:   } else {
    5312:     MatMultAdd(A,x,y,w);
    5313:   }
    5314:   return(0);
    5315: }

    5319: /*@
    5320:    MatInterpolate - y = A*x or A'*x depending on the shape of
    5321:      the matrix

    5323:    Collective on Mat

    5325:    Input Parameters:
    5326: +  mat   - the matrix
    5327: -  x,y - the vectors

    5329:    Level: intermediate

    5331:    Notes:
    5332:     This allows one to use either the restriction or interpolation (its transpose)
    5333:     matrix to do the interpolation

    5335:    Concepts: matrices^interpolation

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

    5339: @*/
    5340: PetscErrorCode PETSCMAT_DLLEXPORT MatInterpolate(Mat A,Vec x,Vec y)
    5341: {
    5343:   PetscInt       M,N;

    5350:   MatPreallocated(A);
    5351:   MatGetSize(A,&M,&N);
    5352:   if (N > M) {
    5353:     MatMultTranspose(A,x,y);
    5354:   } else {
    5355:     MatMult(A,x,y);
    5356:   }
    5357:   return(0);
    5358: }

    5362: /*@
    5363:    MatRestrict - y = A*x or A'*x

    5365:    Collective on Mat

    5367:    Input Parameters:
    5368: +  mat   - the matrix
    5369: -  x,y - the vectors

    5371:    Level: intermediate

    5373:    Notes:
    5374:     This allows one to use either the restriction or interpolation (its transpose)
    5375:     matrix to do the restriction

    5377:    Concepts: matrices^restriction

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

    5381: @*/
    5382: PetscErrorCode PETSCMAT_DLLEXPORT MatRestrict(Mat A,Vec x,Vec y)
    5383: {
    5385:   PetscInt       M,N;

    5392:   MatPreallocated(A);

    5394:   MatGetSize(A,&M,&N);
    5395:   if (N > M) {
    5396:     MatMult(A,x,y);
    5397:   } else {
    5398:     MatMultTranspose(A,x,y);
    5399:   }
    5400:   return(0);
    5401: }

    5405: /*@C
    5406:    MatNullSpaceAttach - attaches a null space to a matrix.
    5407:         This null space will be removed from the resulting vector whenever
    5408:         MatMult() is called

    5410:    Collective on Mat

    5412:    Input Parameters:
    5413: +  mat - the matrix
    5414: -  nullsp - the null space object

    5416:    Level: developer

    5418:    Notes:
    5419:       Overwrites any previous null space that may have been attached

    5421:    Concepts: null space^attaching to matrix

    5423: .seealso: MatCreate(), MatNullSpaceCreate()
    5424: @*/
    5425: PetscErrorCode PETSCMAT_DLLEXPORT MatNullSpaceAttach(Mat mat,MatNullSpace nullsp)
    5426: {

    5433:   MatPreallocated(mat);

    5435:   if (mat->nullsp) {
    5436:     MatNullSpaceDestroy(mat->nullsp);
    5437:   }
    5438:   mat->nullsp = nullsp;
    5439:   PetscObjectReference((PetscObject)nullsp);
    5440:   return(0);
    5441: }

    5445: /*@
    5446:    MatICCFactor - Performs in-place incomplete Cholesky factorization of matrix.

    5448:    Collective on Mat

    5450:    Input Parameters:
    5451: +  mat - the matrix
    5452: .  row - row/column permutation
    5453: .  fill - expected fill factor >= 1.0
    5454: -  level - level of fill, for ICC(k)

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

    5460:    Most users should employ the simplified KSP interface for linear solvers
    5461:    instead of working directly with matrix algebra routines such as this.
    5462:    See, e.g., KSPCreate().

    5464:    Level: developer

    5466:    Concepts: matrices^incomplete Cholesky factorization
    5467:    Concepts: Cholesky factorization

    5469: .seealso: MatICCFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor()
    5470: @*/
    5471: PetscErrorCode PETSCMAT_DLLEXPORT MatICCFactor(Mat mat,IS row,MatFactorInfo* info)
    5472: {

    5480:   if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
    5481:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    5482:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    5483:   if (!mat->ops->iccfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    5484:   MatPreallocated(mat);
    5485:   (*mat->ops->iccfactor)(mat,row,info);
    5486:   PetscObjectStateIncrease((PetscObject)mat);
    5487:   return(0);
    5488: }

    5492: /*@
    5493:    MatSetValuesAdic - Sets values computed with ADIC automatic differentiation into a matrix.

    5495:    Not Collective

    5497:    Input Parameters:
    5498: +  mat - the matrix
    5499: -  v - the values compute with ADIC

    5501:    Level: developer

    5503:    Notes:
    5504:      Must call MatSetColoring() before using this routine. Also this matrix must already
    5505:      have its nonzero pattern determined.

    5507: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
    5508:           MatSetValues(), MatSetColoring(), MatSetValuesAdifor()
    5509: @*/
    5510: PetscErrorCode PETSCMAT_DLLEXPORT MatSetValuesAdic(Mat mat,void *v)
    5511: {

    5519:   if (!mat->assembled) {
    5520:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be already assembled");
    5521:   }
    5522:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
    5523:   if (!mat->ops->setvaluesadic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    5524:   (*mat->ops->setvaluesadic)(mat,v);
    5525:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
    5526:   MatView_Private(mat);
    5527:   PetscObjectStateIncrease((PetscObject)mat);
    5528:   return(0);
    5529: }

    5534: /*@
    5535:    MatSetColoring - Sets a coloring used by calls to MatSetValuesAdic()

    5537:    Not Collective

    5539:    Input Parameters:
    5540: +  mat - the matrix
    5541: -  coloring - the coloring

    5543:    Level: developer

    5545: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
    5546:           MatSetValues(), MatSetValuesAdic()
    5547: @*/
    5548: PetscErrorCode PETSCMAT_DLLEXPORT MatSetColoring(Mat mat,ISColoring coloring)
    5549: {

    5557:   if (!mat->assembled) {
    5558:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be already assembled");
    5559:   }
    5560:   if (!mat->ops->setcoloring) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    5561:   (*mat->ops->setcoloring)(mat,coloring);
    5562:   return(0);
    5563: }

    5567: /*@
    5568:    MatSetValuesAdifor - Sets values computed with automatic differentiation into a matrix.

    5570:    Not Collective

    5572:    Input Parameters:
    5573: +  mat - the matrix
    5574: .  nl - leading dimension of v
    5575: -  v - the values compute with ADIFOR

    5577:    Level: developer

    5579:    Notes:
    5580:      Must call MatSetColoring() before using this routine. Also this matrix must already
    5581:      have its nonzero pattern determined.

    5583: .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(),
    5584:           MatSetValues(), MatSetColoring()
    5585: @*/
    5586: PetscErrorCode PETSCMAT_DLLEXPORT MatSetValuesAdifor(Mat mat,PetscInt nl,void *v)
    5587: {

    5595:   if (!mat->assembled) {
    5596:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be already assembled");
    5597:   }
    5598:   PetscLogEventBegin(MAT_SetValues,mat,0,0,0);
    5599:   if (!mat->ops->setvaluesadifor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    5600:   (*mat->ops->setvaluesadifor)(mat,nl,v);
    5601:   PetscLogEventEnd(MAT_SetValues,mat,0,0,0);
    5602:   PetscObjectStateIncrease((PetscObject)mat);
    5603:   return(0);
    5604: }

    5608: /*@
    5609:    MatDiagonalScaleLocal - Scales columns of a matrix given the scaling values including the
    5610:          ghosted ones.

    5612:    Not Collective

    5614:    Input Parameters:
    5615: +  mat - the matrix
    5616: -  diag = the diagonal values, including ghost ones

    5618:    Level: developer

    5620:    Notes: Works only for MPIAIJ and MPIBAIJ matrices
    5621:
    5622: .seealso: MatDiagonalScale()
    5623: @*/
    5624: PetscErrorCode PETSCMAT_DLLEXPORT MatDiagonalScaleLocal(Mat mat,Vec diag)
    5625: {
    5627:   PetscMPIInt    size;

    5634:   if (!mat->assembled) {
    5635:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be already assembled");
    5636:   }
    5637:   PetscLogEventBegin(MAT_Scale,mat,0,0,0);
    5638:   MPI_Comm_size(mat->comm,&size);
    5639:   if (size == 1) {
    5640:     PetscInt n,m;
    5641:     VecGetSize(diag,&n);
    5642:     MatGetSize(mat,0,&m);
    5643:     if (m == n) {
    5644:       MatDiagonalScale(mat,0,diag);
    5645:     } else {
    5646:       SETERRQ(PETSC_ERR_SUP,"Only supported for sequential matrices when no ghost points/periodic conditions");
    5647:     }
    5648:   } else {
    5649:     PetscErrorCode (*f)(Mat,Vec);
    5650:     PetscObjectQueryFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",(void (**)(void))&f);
    5651:     if (f) {
    5652:       (*f)(mat,diag);
    5653:     } else {
    5654:       SETERRQ(PETSC_ERR_SUP,"Only supported for MPIAIJ and MPIBAIJ parallel matrices");
    5655:     }
    5656:   }
    5657:   PetscLogEventEnd(MAT_Scale,mat,0,0,0);
    5658:   PetscObjectStateIncrease((PetscObject)mat);
    5659:   return(0);
    5660: }

    5664: /*@
    5665:    MatGetInertia - Gets the inertia from a factored matrix

    5667:    Collective on Mat

    5669:    Input Parameter:
    5670: .  mat - the matrix

    5672:    Output Parameters:
    5673: +   nneg - number of negative eigenvalues
    5674: .   nzero - number of zero eigenvalues
    5675: -   npos - number of positive eigenvalues

    5677:    Level: advanced

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

    5682: @*/
    5683: PetscErrorCode PETSCMAT_DLLEXPORT MatGetInertia(Mat mat,PetscInt *nneg,PetscInt *nzero,PetscInt *npos)
    5684: {

    5690:   if (!mat->factor)    SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
    5691:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Numeric factor mat is not assembled");
    5692:   if (!mat->ops->getinertia) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    5693:   (*mat->ops->getinertia)(mat,nneg,nzero,npos);
    5694:   return(0);
    5695: }

    5697: /* ----------------------------------------------------------------*/
    5700: /*@
    5701:    MatSolves - Solves A x = b, given a factored matrix, for a collection of vectors

    5703:    Collective on Mat and Vecs

    5705:    Input Parameters:
    5706: +  mat - the factored matrix
    5707: -  b - the right-hand-side vectors

    5709:    Output Parameter:
    5710: .  x - the result vectors

    5712:    Notes:
    5713:    The vectors b and x cannot be the same.  I.e., one cannot
    5714:    call MatSolves(A,x,x).

    5716:    Notes:
    5717:    Most users should employ the simplified KSP interface for linear solvers
    5718:    instead of working directly with matrix algebra routines such as this.
    5719:    See, e.g., KSPCreate().

    5721:    Level: developer

    5723:    Concepts: matrices^triangular solves

    5725: .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd(), MatSolve()
    5726: @*/
    5727: PetscErrorCode PETSCMAT_DLLEXPORT MatSolves(Mat mat,Vecs b,Vecs x)
    5728: {

    5734:   if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors");
    5735:   if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix");
    5736:   if (!mat->M && !mat->N) return(0);

    5738:   if (!mat->ops->solves) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name);
    5739:   MatPreallocated(mat);
    5740:   PetscLogEventBegin(MAT_Solves,mat,0,0,0);
    5741:   (*mat->ops->solves)(mat,b,x);
    5742:   PetscLogEventEnd(MAT_Solves,mat,0,0,0);
    5743:   return(0);
    5744: }

    5748: /*@
    5749:    MatIsSymmetric - Test whether a matrix is symmetric

    5751:    Collective on Mat

    5753:    Input Parameter:
    5754: +  A - the matrix to test
    5755: -  tol - difference between value and its transpose less than this amount counts as equal (use 0.0 for exact transpose)

    5757:    Output Parameters:
    5758: .  flg - the result

    5760:    Level: intermediate

    5762:    Concepts: matrix^symmetry

    5764: .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsStructurallySymmetric(), MatSetOption(), MatIsSymmetricKnown()
    5765: @*/
    5766: PetscErrorCode PETSCMAT_DLLEXPORT MatIsSymmetric(Mat A,PetscReal tol,PetscTruth *flg)
    5767: {

    5773:   if (!A->symmetric_set) {
    5774:     if (!A->ops->issymmetric) {
    5775:       MatType mattype;
    5776:       MatGetType(A,&mattype);
    5777:       SETERRQ1(PETSC_ERR_SUP,"Matrix of type <%s> does not support checking for symmetric",mattype);
    5778:     }
    5779:     (*A->ops->issymmetric)(A,tol,&A->symmetric);
    5780:     A->symmetric_set = PETSC_TRUE;
    5781:     if (A->symmetric) {
    5782:       A->structurally_symmetric_set = PETSC_TRUE;
    5783:       A->structurally_symmetric     = PETSC_TRUE;
    5784:     }
    5785:   }
    5786:   *flg = A->symmetric;
    5787:   return(0);
    5788: }

    5792: /*@
    5793:    MatIsSymmetricKnown - Checks the flag on the matrix to see if it is symmetric.

    5795:    Collective on Mat

    5797:    Input Parameter:
    5798: .  A - the matrix to check

    5800:    Output Parameters:
    5801: +  set - if the symmetric flag is set (this tells you if the next flag is valid)
    5802: -  flg - the result

    5804:    Level: advanced

    5806:    Concepts: matrix^symmetry

    5808:    Note: Does not check the matrix values directly, so this may return unknown (set = PETSC_FALSE). Use MatIsSymmetric()
    5809:          if you want it explicitly checked

    5811: .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsStructurallySymmetric(), MatSetOption(), MatIsSymmetric()
    5812: @*/
    5813: PetscErrorCode PETSCMAT_DLLEXPORT MatIsSymmetricKnown(Mat A,PetscTruth *set,PetscTruth *flg)
    5814: {
    5819:   if (A->symmetric_set) {
    5820:     *set = PETSC_TRUE;
    5821:     *flg = A->symmetric;
    5822:   } else {
    5823:     *set = PETSC_FALSE;
    5824:   }
    5825:   return(0);
    5826: }

    5830: /*@
    5831:    MatIsHermitianKnown - Checks the flag on the matrix to see if it is hermitian.

    5833:    Collective on Mat

    5835:    Input Parameter:
    5836: .  A - the matrix to check

    5838:    Output Parameters:
    5839: +  set - if the hermitian flag is set (this tells you if the next flag is valid)
    5840: -  flg - the result

    5842:    Level: advanced

    5844:    Concepts: matrix^symmetry

    5846:    Note: Does not check the matrix values directly, so this may return unknown (set = PETSC_FALSE). Use MatIsHermitian()
    5847:          if you want it explicitly checked

    5849: .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsStructurallySymmetric(), MatSetOption(), MatIsSymmetric()
    5850: @*/
    5851: PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianKnown(Mat A,PetscTruth *set,PetscTruth *flg)
    5852: {
    5857:   if (A->hermitian_set) {
    5858:     *set = PETSC_TRUE;
    5859:     *flg = A->hermitian;
    5860:   } else {
    5861:     *set = PETSC_FALSE;
    5862:   }
    5863:   return(0);
    5864: }

    5868: /*@
    5869:    MatIsStructurallySymmetric - Test whether a matrix is structurally symmetric

    5871:    Collective on Mat

    5873:    Input Parameter:
    5874: .  A - the matrix to test

    5876:    Output Parameters:
    5877: .  flg - the result

    5879:    Level: intermediate

    5881:    Concepts: matrix^symmetry

    5883: .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsSymmetric(), MatSetOption()
    5884: @*/
    5885: PetscErrorCode PETSCMAT_DLLEXPORT MatIsStructurallySymmetric(Mat A,PetscTruth *flg)
    5886: {

    5892:   if (!A->structurally_symmetric_set) {
    5893:     if (!A->ops->isstructurallysymmetric) SETERRQ(PETSC_ERR_SUP,"Matrix does not support checking for structural symmetric");
    5894:     (*A->ops->isstructurallysymmetric)(A,&A->structurally_symmetric);
    5895:     A->structurally_symmetric_set = PETSC_TRUE;
    5896:   }
    5897:   *flg = A->structurally_symmetric;
    5898:   return(0);
    5899: }

    5903: /*@
    5904:    MatIsHermitian - Test whether a matrix is Hermitian, i.e. it is the complex conjugate of its transpose.

    5906:    Collective on Mat

    5908:    Input Parameter:
    5909: .  A - the matrix to test

    5911:    Output Parameters:
    5912: .  flg - the result

    5914:    Level: intermediate

    5916:    Concepts: matrix^symmetry

    5918: .seealso: MatTranspose(), MatIsTranspose(), MatIsSymmetric(), MatIsStructurallySymmetric(), MatSetOption()
    5919: @*/
    5920: PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitian(Mat A,PetscTruth *flg)
    5921: {

    5927:   if (!A->hermitian_set) {
    5928:     if (!A->ops->ishermitian) SETERRQ(PETSC_ERR_SUP,"Matrix does not support checking for being Hermitian");
    5929:     (*A->ops->ishermitian)(A,&A->hermitian);
    5930:     A->hermitian_set = PETSC_TRUE;
    5931:     if (A->hermitian) {
    5932:       A->structurally_symmetric_set = PETSC_TRUE;
    5933:       A->structurally_symmetric     = PETSC_TRUE;
    5934:     }
    5935:   }
    5936:   *flg = A->hermitian;
    5937:   return(0);
    5938: }

    5943: /*@
    5944:    MatStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
    5945:        to be communicated to other processors during the MatAssemblyBegin/End() process

    5947:     Not collective

    5949:    Input Parameter:
    5950: .   vec - the vector

    5952:    Output Parameters:
    5953: +   nstash   - the size of the stash
    5954: .   reallocs - the number of additional mallocs incurred.
    5955: .   bnstash   - the size of the block stash
    5956: -   breallocs - the number of additional mallocs incurred.in the block stash
    5957:
    5958:    Level: advanced

    5960: .seealso: MatAssemblyBegin(), MatAssemblyEnd(), Mat, MatStashSetInitialSize()
    5961:
    5962: @*/
    5963: PetscErrorCode PETSCMAT_DLLEXPORT MatStashGetInfo(Mat mat,PetscInt *nstash,PetscInt *reallocs,PetscInt *bnstash,PetscInt *brealloc)
    5964: {
    5967:   MatStashGetInfo_Private(&mat->stash,nstash,reallocs);
    5968:   MatStashGetInfo_Private(&mat->bstash,nstash,reallocs);
    5969:   return(0);
    5970: }

    5974: /*@
    5975:    MatGetVecs - Get vector(s) compatible with the matrix, i.e. with the same
    5976:      parallel layout
    5977:
    5978:    Collective on Mat

    5980:    Input Parameter:
    5981: .  mat - the matrix

    5983:    Output Parameter:
    5984: +   right - (optional) vector that the matrix can be multiplied against
    5985: -   left - (optional) vector that the matrix vector product can be stored in

    5987:   Level: advanced

    5989: .seealso: MatCreate()
    5990: @*/
    5991: PetscErrorCode PETSCMAT_DLLEXPORT MatGetVecs(Mat mat,Vec *right,Vec *left)
    5992: {

    5998:   MatPreallocated(mat);
    5999:   if (mat->ops->getvecs) {
    6000:     (*mat->ops->getvecs)(mat,right,left);
    6001:   } else {
    6002:     PetscMPIInt size;
    6003:     MPI_Comm_size(mat->comm, &size);
    6004:     if (right) {
    6005:       VecCreate(mat->comm,right);
    6006:       VecSetSizes(*right,mat->n,PETSC_DETERMINE);
    6007:       if (size > 1) {VecSetType(*right,VECMPI);}
    6008:       else {VecSetType(*right,VECSEQ);}
    6009:     }
    6010:     if (left) {
    6011:       VecCreate(mat->comm,left);
    6012:       VecSetSizes(*left,mat->m,PETSC_DETERMINE);
    6013:       if (size > 1) {VecSetType(*left,VECMPI);}
    6014:       else {VecSetType(*left,VECSEQ);}
    6015:     }
    6016:   }
    6017:   if (right) {VecSetBlockSize(*right,mat->bs);}
    6018:   if (left) {VecSetBlockSize(*left,mat->bs);}
    6019:   return(0);
    6020: }

    6024: /*@
    6025:    MatFactorInfoInitialize - Initializes a MatFactorInfo data structure
    6026:      with default values.

    6028:    Not Collective

    6030:    Input Parameters:
    6031: .    info - the MatFactorInfo data structure

    6034:    Notes: The solvers are generally used through the KSP and PC objects, for example
    6035:           PCLU, PCILU, PCCHOLESKY, PCICC

    6037:    Level: developer

    6039: .seealso: MatFactorInfo
    6040: @*/

    6042: PetscErrorCode PETSCMAT_DLLEXPORT MatFactorInfoInitialize(MatFactorInfo *info)
    6043: {

    6047:   PetscMemzero(info,sizeof(MatFactorInfo));
    6048:   return(0);
    6049: }

    6053: /*@
    6054:    MatPtAP - Creates the matrix projection C = P^T * A * P

    6056:    Collective on Mat

    6058:    Input Parameters:
    6059: +  A - the matrix
    6060: .  P - the projection matrix
    6061: .  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
    6062: -  fill - expected fill as ratio of nnz(C)/(nnz(A) + nnz(P))

    6064:    Output Parameters:
    6065: .  C - the product matrix

    6067:    Notes:
    6068:    C will be created and must be destroyed by the user with MatDestroy().

    6070:    This routine is currently only implemented for pairs of AIJ matrices and classes
    6071:    which inherit from AIJ.

    6073:    Level: intermediate

    6075: .seealso: MatPtAPSymbolic(),MatPtAPNumeric(),MatMatMult()
    6076: @*/
    6077: PetscErrorCode PETSCMAT_DLLEXPORT MatPtAP(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
    6078: {

    6084:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6085:   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    6088:   MatPreallocated(P);
    6089:   if (!P->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6090:   if (P->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    6092:   if (P->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",P->M,A->N);
    6093:   if (fill <=0.0) SETERRQ1(PETSC_ERR_ARG_SIZ,"fill=%g must be > 0.0",fill);
    6094:   MatPreallocated(A);

    6096:   PetscLogEventBegin(MAT_PtAP,A,P,0,0);
    6097:   (*A->ops->ptap)(A,P,scall,fill,C);
    6098:   PetscLogEventEnd(MAT_PtAP,A,P,0,0);

    6100:   return(0);
    6101: }

    6105: /*@
    6106:    MatPtAPNumeric - Computes the matrix projection C = P^T * A * P

    6108:    Collective on Mat

    6110:    Input Parameters:
    6111: +  A - the matrix
    6112: -  P - the projection matrix

    6114:    Output Parameters:
    6115: .  C - the product matrix

    6117:    Notes:
    6118:    C must have been created by calling MatPtAPSymbolic and must be destroyed by
    6119:    the user using MatDeatroy().

    6121:    This routine is currently only implemented for pairs of AIJ matrices and classes
    6122:    which inherit from AIJ.  C will be of type MATAIJ.

    6124:    Level: intermediate

    6126: .seealso: MatPtAP(),MatPtAPSymbolic(),MatMatMultNumeric()
    6127: @*/
    6128: PetscErrorCode PETSCMAT_DLLEXPORT MatPtAPNumeric(Mat A,Mat P,Mat C)
    6129: {

    6135:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6136:   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    6139:   MatPreallocated(P);
    6140:   if (!P->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6141:   if (P->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    6144:   MatPreallocated(C);
    6145:   if (C->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    6146:   if (P->N!=C->M) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",P->N,C->M);
    6147:   if (P->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",P->M,A->N);
    6148:   if (A->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix 'A' must be square, %D != %D",A->M,A->N);
    6149:   if (P->N!=C->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",P->N,C->N);
    6150:   MatPreallocated(A);

    6152:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
    6153:   (*A->ops->ptapnumeric)(A,P,C);
    6154:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
    6155:   return(0);
    6156: }

    6160: /*@
    6161:    MatPtAPSymbolic - Creates the (i,j) structure of the matrix projection C = P^T * A * P

    6163:    Collective on Mat

    6165:    Input Parameters:
    6166: +  A - the matrix
    6167: -  P - the projection matrix

    6169:    Output Parameters:
    6170: .  C - the (i,j) structure of the product matrix

    6172:    Notes:
    6173:    C will be created and must be destroyed by the user with MatDestroy().

    6175:    This routine is currently only implemented for pairs of SeqAIJ matrices and classes
    6176:    which inherit from SeqAIJ.  C will be of type MATSEQAIJ.  The product is computed using
    6177:    this (i,j) structure by calling MatPtAPNumeric().

    6179:    Level: intermediate

    6181: .seealso: MatPtAP(),MatPtAPNumeric(),MatMatMultSymbolic()
    6182: @*/
    6183: PetscErrorCode PETSCMAT_DLLEXPORT MatPtAPSymbolic(Mat A,Mat P,PetscReal fill,Mat *C)
    6184: {

    6190:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6191:   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    6194:   MatPreallocated(P);
    6195:   if (!P->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6196:   if (P->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

    6199:   if (P->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",P->M,A->N);
    6200:   if (A->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix 'A' must be square, %D != %D",A->M,A->N);
    6201:   MatPreallocated(A);
    6202:   PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
    6203:   (*A->ops->ptapsymbolic)(A,P,fill,C);
    6204:   PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);

    6206:   MatSetBlockSize(*C,A->bs);

    6208:   return(0);
    6209: }

    6213: /*@
    6214:    MatMatMult - Performs Matrix-Matrix Multiplication C=A*B.

    6216:    Collective on Mat

    6218:    Input Parameters:
    6219: +  A - the left matrix
    6220: .  B - the right matrix
    6221: .  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
    6222: -  fill - expected fill as ratio of nnz(C)/(nnz(A) + nnz(B))

    6224:    Output Parameters:
    6225: .  C - the product matrix

    6227:    Notes:
    6228:    C will be created and must be destroyed by the user with MatDestroy().

    6230:    This routine is currently only implemented for pairs of AIJ matrices and classes
    6231:    which inherit from AIJ.  C will be of type MATAIJ.

    6233:    Level: intermediate

    6235: .seealso: MatMatMultSymbolic(),MatMatMultNumeric()
    6236: @*/
    6237: PetscErrorCode PETSCMAT_DLLEXPORT MatMatMult(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
    6238: {
    6240:   PetscErrorCode (*fA)(Mat,Mat,MatReuse,PetscReal,Mat*);
    6241:   PetscErrorCode (*fB)(Mat,Mat,MatReuse,PetscReal,Mat*);

    6246:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6247:   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    6250:   MatPreallocated(B);
    6251:   if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6252:   if (B->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    6254:   if (B->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->M,A->N);
    6255:   if (fill <=0.0) SETERRQ1(PETSC_ERR_ARG_SIZ,"fill=%g must be > 0.0",fill);
    6256:   MatPreallocated(A);

    6258:   /* For now, we do not dispatch based on the type of A and B */
    6259:   /* When implementations like _SeqAIJ_MAIJ exist, attack the multiple dispatch problem. */
    6260:   fA = A->ops->matmult;
    6261:   if (!fA) SETERRQ1(PETSC_ERR_SUP,"MatMatMult not supported for A of type %s",A->type_name);
    6262:   fB = B->ops->matmult;
    6263:   if (!fB) SETERRQ1(PETSC_ERR_SUP,"MatMatMult not supported for B of type %s",B->type_name);
    6264:   if (fB!=fA) SETERRQ2(PETSC_ERR_ARG_INCOMP,"MatMatMult requires A, %s, to be compatible with B, %s",A->type_name,B->type_name);

    6266:   PetscLogEventBegin(MAT_MatMult,A,B,0,0);
    6267:   (*A->ops->matmult)(A,B,scall,fill,C);
    6268:   PetscLogEventEnd(MAT_MatMult,A,B,0,0);
    6269:
    6270:   return(0);
    6271: }

    6275: /*@
    6276:    MatMatMultSymbolic - Performs construction, preallocation, and computes the ij structure
    6277:    of the matrix-matrix product C=A*B.  Call this routine before calling MatMatMultNumeric().

    6279:    Collective on Mat

    6281:    Input Parameters:
    6282: +  A - the left matrix
    6283: .  B - the right matrix
    6284: -  fill - expected fill as ratio of nnz(C)/(nnz(A) + nnz(B))

    6286:    Output Parameters:
    6287: .  C - the matrix containing the ij structure of product matrix

    6289:    Notes:
    6290:    C will be created as a MATSEQAIJ matrix and must be destroyed by the user with MatDestroy().

    6292:    This routine is currently only implemented for SeqAIJ matrices and classes which inherit from SeqAIJ.

    6294:    Level: intermediate

    6296: .seealso: MatMatMult(),MatMatMultNumeric()
    6297: @*/
    6298: PetscErrorCode PETSCMAT_DLLEXPORT MatMatMultSymbolic(Mat A,Mat B,PetscReal fill,Mat *C)
    6299: {
    6301:   PetscErrorCode (*Asymbolic)(Mat,Mat,PetscReal,Mat *);
    6302:   PetscErrorCode (*Bsymbolic)(Mat,Mat,PetscReal,Mat *);

    6307:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6308:   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

    6312:   MatPreallocated(B);
    6313:   if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6314:   if (B->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

    6317:   if (B->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->M,A->N);
    6318:   if (fill <=0.0) SETERRQ1(PETSC_ERR_ARG_SIZ,"fill=%g must be > 0.0",fill);
    6319:   MatPreallocated(A);

    6321:   /* For now, we do not dispatch based on the type of A and P */
    6322:   /* When implementations like _SeqAIJ_MAIJ exist, attack the multiple dispatch problem. */
    6323:   Asymbolic = A->ops->matmultsymbolic;
    6324:   if (!Asymbolic) SETERRQ1(PETSC_ERR_SUP,"C=A*B not implemented for A of type %s",A->type_name);
    6325:   Bsymbolic = B->ops->matmultsymbolic;
    6326:   if (!Bsymbolic) SETERRQ1(PETSC_ERR_SUP,"C=A*B not implemented for B of type %s",B->type_name);
    6327:   if (Bsymbolic!=Asymbolic) SETERRQ2(PETSC_ERR_ARG_INCOMP,"MatMatMultSymbolic requires A, %s, to be compatible with B, %s",A->type_name,B->type_name);

    6329:   PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
    6330:   (*Asymbolic)(A,B,fill,C);
    6331:   PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);

    6333:   return(0);
    6334: }

    6338: /*@
    6339:    MatMatMultNumeric - Performs the numeric matrix-matrix product.
    6340:    Call this routine after first calling MatMatMultSymbolic().

    6342:    Collective on Mat

    6344:    Input Parameters:
    6345: +  A - the left matrix
    6346: -  B - the right matrix

    6348:    Output Parameters:
    6349: .  C - the product matrix, whose ij structure was defined from MatMatMultSymbolic().

    6351:    Notes:
    6352:    C must have been created with MatMatMultSymbolic.

    6354:    This routine is currently only implemented for SeqAIJ type matrices.

    6356:    Level: intermediate

    6358: .seealso: MatMatMult(),MatMatMultSymbolic()
    6359: @*/
    6360: PetscErrorCode PETSCMAT_DLLEXPORT MatMatMultNumeric(Mat A,Mat B,Mat C)
    6361: {
    6363:   PetscErrorCode (*Anumeric)(Mat,Mat,Mat);
    6364:   PetscErrorCode (*Bnumeric)(Mat,Mat,Mat);

    6370:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6371:   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

    6375:   MatPreallocated(B);
    6376:   if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6377:   if (B->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

    6381:   MatPreallocated(C);
    6382:   if (!C->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6383:   if (C->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

    6385:   if (B->N!=C->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->N,C->N);
    6386:   if (B->M!=A->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->M,A->N);
    6387:   if (A->M!=C->M) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",A->M,C->M);
    6388:   MatPreallocated(A);

    6390:   /* For now, we do not dispatch based on the type of A and B */
    6391:   /* When implementations like _SeqAIJ_MAIJ exist, attack the multiple dispatch problem. */
    6392:   Anumeric = A->ops->matmultnumeric;
    6393:   if (!Anumeric) SETERRQ1(PETSC_ERR_SUP,"MatMatMultNumeric not supported for A of type %s",A->type_name);
    6394:   Bnumeric = B->ops->matmultnumeric;
    6395:   if (!Bnumeric) SETERRQ1(PETSC_ERR_SUP,"MatMatMultNumeric not supported for B of type %s",B->type_name);
    6396:   if (Bnumeric!=Anumeric) SETERRQ2(PETSC_ERR_ARG_INCOMP,"MatMatMultNumeric requires A, %s, to be compatible with B, %s",A->type_name,B->type_name);

    6398:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
    6399:   (*Anumeric)(A,B,C);
    6400:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);

    6402:   return(0);
    6403: }

    6407: /*@
    6408:    MatMatMultTranspose - Performs Matrix-Matrix Multiplication C=A^T*B.

    6410:    Collective on Mat

    6412:    Input Parameters:
    6413: +  A - the left matrix
    6414: .  B - the right matrix
    6415: .  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
    6416: -  fill - expected fill as ratio of nnz(C)/(nnz(A) + nnz(B))

    6418:    Output Parameters:
    6419: .  C - the product matrix

    6421:    Notes:
    6422:    C will be created and must be destroyed by the user with MatDestroy().

    6424:    This routine is currently only implemented for pairs of SeqAIJ matrices and classes
    6425:    which inherit from SeqAIJ.  C will be of type MATSEQAIJ.

    6427:    Level: intermediate

    6429: .seealso: MatMatMultTransposeSymbolic(),MatMatMultTransposeNumeric()
    6430: @*/
    6431: PetscErrorCode PETSCMAT_DLLEXPORT MatMatMultTranspose(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
    6432: {
    6434:   PetscErrorCode (*fA)(Mat,Mat,MatReuse,PetscReal,Mat*);
    6435:   PetscErrorCode (*fB)(Mat,Mat,MatReuse,PetscReal,Mat*);

    6440:   if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6441:   if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    6444:   MatPreallocated(B);
    6445:   if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    6446:   if (B->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
    6448:   if (B->M!=A->M) SETERRQ2(PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->M,A->M);
    6449:   if (fill <=0.0) SETERRQ1(PETSC_ERR_ARG_SIZ,"fill=%g must be > 0.0",fill);
    6450:   MatPreallocated(A);

    6452:   fA = A->ops->matmulttranspose;
    6453:   if (!fA) SETERRQ1(PETSC_ERR_SUP,"MatMatMultTranspose not supported for A of type %s",A->type_name);
    6454:   fB = B->ops->matmulttranspose;
    6455:   if (!fB) SETERRQ1(PETSC_ERR_SUP,"MatMatMultTranspose not supported for B of type %s",B->type_name);
    6456:   if (fB!=fA) SETERRQ2(PETSC_ERR_ARG_INCOMP,"MatMatMultTranspose requires A, %s, to be compatible with B, %s",A->type_name,B->type_name);

    6458:   PetscLogEventBegin(MAT_MatMultTranspose,A,B,0,0);
    6459:   (*A->ops->matmulttranspose)(A,B,scall,fill,C);
    6460:   PetscLogEventEnd(MAT_MatMultTranspose,A,B,0,0);
    6461:
    6462:   return(0);
    6463: }