Actual source code: matreg.c

      1: /*
      2:      Mechanism for register PETSc matrix types
      3: */
      4: #include <petsc/private/matimpl.h>

      6: PetscBool MatRegisterAllCalled = PETSC_FALSE;

      8: /*
      9:    Contains the list of registered Mat routines
     10: */
     11: PetscFunctionList MatList = NULL;

     13: /* MatGetRootType_Private - Gets the root type of the input matrix's type (e.g., MATAIJ for MATSEQAIJ)

     15:    Not Collective

     17:    Input Parameter:
     18: .  mat      - the input matrix, could be sequential or MPI

     20:    Output Parameter:
     21: .  rootType  - the root matrix type

     23:    Level: developer

     25: .seealso: `MatGetType()`, `MatSetType()`, `MatType`, `Mat`
     26: */
     27: PetscErrorCode MatGetRootType_Private(Mat mat, MatType *rootType)
     28: {
     29:   PetscBool   found = PETSC_FALSE;
     30:   MatRootName names = MatRootNameList;
     31:   MatType     inType;

     33:   PetscFunctionBegin;
     35:   PetscCall(MatGetType(mat, &inType));
     36:   while (names) {
     37:     PetscCall(PetscStrcmp(inType, names->mname, &found));
     38:     if (!found) PetscCall(PetscStrcmp(inType, names->sname, &found));
     39:     if (found) {
     40:       found     = PETSC_TRUE;
     41:       *rootType = names->rname;
     42:       break;
     43:     }
     44:     names = names->next;
     45:   }
     46:   if (!found) *rootType = inType;
     47:   PetscFunctionReturn(PETSC_SUCCESS);
     48: }

     50: /* MatGetMPIMatType_Private - Gets the MPI type corresponding to the input matrix's type (e.g., MATMPIAIJ for MATSEQAIJ)

     52:    Not Collective

     54:    Input Parameter:
     55: .  mat      - the input matrix, could be sequential or MPI

     57:    Output Parameter:
     58: .  MPIType  - the parallel (MPI) matrix type

     60:    Level: developer

     62: .seealso: `MatGetType()`, `MatSetType()`, `MatType`, `Mat`
     63: */
     64: PetscErrorCode MatGetMPIMatType_Private(Mat mat, MatType *MPIType)
     65: {
     66:   PetscBool   found = PETSC_FALSE;
     67:   MatRootName names = MatRootNameList;
     68:   MatType     inType;

     70:   PetscFunctionBegin;
     72:   PetscCall(MatGetType(mat, &inType));
     73:   while (names) {
     74:     PetscCall(PetscStrcmp(inType, names->sname, &found));
     75:     if (!found) PetscCall(PetscStrcmp(inType, names->mname, &found));
     76:     if (!found) PetscCall(PetscStrcmp(inType, names->rname, &found));
     77:     if (found) {
     78:       found    = PETSC_TRUE;
     79:       *MPIType = names->mname;
     80:       break;
     81:     }
     82:     names = names->next;
     83:   }
     84:   PetscCheck(found, PETSC_COMM_SELF, PETSC_ERR_SUP, "No corresponding parallel (MPI) type for this matrix");
     85:   PetscFunctionReturn(PETSC_SUCCESS);
     86: }

     88: /*@
     89:   MatSetType - Builds matrix object for a particular matrix type

     91:   Collective

     93:   Input Parameters:
     94: + mat    - the matrix object
     95: - matype - matrix type

     97:   Options Database Key:
     98: . -mat_type type - Sets the type; see `MatType`

    100:   Level: intermediate

    102:   Note:
    103:   See `MatType` for possible values

    105: .seealso: [](ch_matrices), `Mat`, `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`
    106: @*/
    107: PetscErrorCode MatSetType(Mat mat, MatType matype)
    108: {
    109:   PetscBool   sametype, found, subclass = PETSC_FALSE, matMPI = PETSC_FALSE, requestSeq = PETSC_FALSE;
    110:   MatRootName names = MatRootNameList;
    111:   PetscErrorCode (*r)(Mat);

    113:   PetscFunctionBegin;

    116:   /* make this special case fast */
    117:   PetscCall(PetscObjectTypeCompare((PetscObject)mat, matype, &sametype));
    118:   if (sametype) PetscFunctionReturn(PETSC_SUCCESS);

    120:   /* suppose with one MPI process, one created an MPIAIJ (mpiaij) matrix with MatCreateMPIAIJWithArrays(), and later tried
    121:      to change its type via '-mat_type aijcusparse'. Even there is only one MPI rank, we need to adapt matype to
    122:      'mpiaijcusparse' so that it will be treated as a subclass of MPIAIJ and proper MatCovert() will be called.
    123:   */
    124:   if (((PetscObject)mat)->type_name) PetscCall(PetscStrbeginswith(((PetscObject)mat)->type_name, "mpi", &matMPI)); /* mat claims itself is an 'mpi' matrix */
    125:   if (matype) PetscCall(PetscStrbeginswith(matype, "seq", &requestSeq));                                           /* user is requesting a 'seq' matrix */
    126:   PetscCheck(!(matMPI && requestSeq), PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Changing an MPI matrix (%s) to a sequential one (%s) is not allowed. Please remove the 'seq' prefix from your matrix type.", ((PetscObject)mat)->type_name, matype);

    128:   while (names) {
    129:     PetscCall(PetscStrcmp(matype, names->rname, &found));
    130:     if (found) {
    131:       PetscMPIInt size;
    132:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size));
    133:       if (size == 1 && !matMPI) matype = names->sname; /* try to align the requested type (matype) with the existing type per seq/mpi */
    134:       else matype = names->mname;
    135:       break;
    136:     }
    137:     names = names->next;
    138:   }

    140:   PetscCall(PetscObjectTypeCompare((PetscObject)mat, matype, &sametype));
    141:   if (sametype) PetscFunctionReturn(PETSC_SUCCESS);

    143:   PetscCall(PetscFunctionListFind(MatList, matype, &r));
    144:   PetscCheck(r, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown Mat type given: %s", matype);

    146:   if (mat->assembled && ((PetscObject)mat)->type_name) PetscCall(PetscStrbeginswith(matype, ((PetscObject)mat)->type_name, &subclass));
    147:   if (subclass) { /* mat is a subclass of the requested 'matype'? */
    148:     PetscCall(MatConvert(mat, matype, MAT_INPLACE_MATRIX, &mat));
    149:     PetscFunctionReturn(PETSC_SUCCESS);
    150:   }
    151:   if (names && mat->assembled) {
    152:     PetscCall(PetscStrbeginswith(names->rname, "sell", &sametype));
    153:     if (sametype) {                                                  /* mattype is MATSELL or its subclass */
    154:       PetscCall(MatConvert(mat, MATSELL, MAT_INPLACE_MATRIX, &mat)); /* convert to matsell first */
    155:       PetscCall(MatConvert(mat, matype, MAT_INPLACE_MATRIX, &mat));
    156:       PetscFunctionReturn(PETSC_SUCCESS);
    157:     }
    158:   }
    159:   PetscTryTypeMethod(mat, destroy);
    160:   mat->ops->destroy = NULL;

    162:   /* should these null spaces be removed? */
    163:   PetscCall(MatNullSpaceDestroy(&mat->nullsp));
    164:   PetscCall(MatNullSpaceDestroy(&mat->nearnullsp));

    166:   PetscCall(PetscMemzero(mat->ops, sizeof(struct _MatOps)));
    167:   mat->preallocated  = PETSC_FALSE;
    168:   mat->assembled     = PETSC_FALSE;
    169:   mat->was_assembled = PETSC_FALSE;

    171:   /*
    172:    Increment, rather than reset these: the object is logically the same, so its logging and
    173:    state is inherited.  Furthermore, resetting makes it possible for the same state to be
    174:    obtained with a different structure, confusing the PC.
    175:   */
    176:   mat->nonzerostate++;
    177:   PetscCall(PetscObjectStateIncrease((PetscObject)mat));

    179:   /* create the new data structure */
    180:   PetscCall((*r)(mat));
    181:   PetscFunctionReturn(PETSC_SUCCESS);
    182: }

    184: /*@
    185:   MatGetType - Gets the matrix type as a string from the matrix object.

    187:   Not Collective

    189:   Input Parameter:
    190: . mat - the matrix

    192:   Output Parameter:
    193: . type - name of matrix type

    195:   Level: intermediate

    197:   Note:
    198:   `type` should not be retained for later use as it will be an invalid pointer if the `MatType` of `mat` is changed.

    200: .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetType()`, `PetscObjectTypeCompare()`, `PetscObjectTypeCompareAny()`
    201: @*/
    202: PetscErrorCode MatGetType(Mat mat, MatType *type)
    203: {
    204:   PetscFunctionBegin;
    206:   PetscAssertPointer(type, 2);
    207:   *type = ((PetscObject)mat)->type_name;
    208:   PetscFunctionReturn(PETSC_SUCCESS);
    209: }

    211: /*@
    212:   MatGetVecType - Gets the vector type the matrix will return with `MatCreateVecs()`

    214:   Not Collective

    216:   Input Parameter:
    217: . mat - the matrix

    219:   Output Parameter:
    220: . vtype - name of vector type

    222:   Level: intermediate

    224: .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetVecType()`, `VecType`
    225: @*/
    226: PetscErrorCode MatGetVecType(Mat mat, VecType *vtype)
    227: {
    228:   PetscFunctionBegin;
    230:   PetscAssertPointer(vtype, 2);
    231:   *vtype = mat->defaultvectype;
    232:   PetscFunctionReturn(PETSC_SUCCESS);
    233: }

    235: /*@
    236:   MatSetVecType - Set the vector type the matrix will return with `MatCreateVecs()`

    238:   Collective

    240:   Input Parameters:
    241: + mat   - the matrix object
    242: - vtype - vector type

    244:   Options Database Key:
    245: . -mat_vec_type vtype - set the `VecType` of the created vectors during `MatSetFromOptions()`

    247:   Level: advanced

    249:   Note:
    250:   This is rarely needed in practice since each matrix object internally sets the proper vector type.

    252: .seealso: [](ch_matrices), `Mat`, `VecType`, `VecSetType()`, `MatGetVecType()`, `MatCreateVecs()`
    253: @*/
    254: PetscErrorCode MatSetVecType(Mat mat, VecType vtype)
    255: {
    256:   PetscFunctionBegin;
    258:   if (vtype != mat->defaultvectype) {
    259:     PetscCall(PetscFree(mat->defaultvectype));
    260:     PetscCall(PetscStrallocpy(vtype, &mat->defaultvectype));
    261:   }
    262:   PetscFunctionReturn(PETSC_SUCCESS);
    263: }

    265: /*@C
    266:   MatRegister - Adds a new matrix type implementation that is usable as a `Mat` in PETSc

    268:   Not Collective, No Fortran Support

    270:   Input Parameters:
    271: + sname    - name of a new user-defined matrix type
    272: - function - routine to create method context

    274:   Level: advanced

    276:   Note:
    277:   `MatRegister()` may be called multiple times to add several user-defined solvers.

    279:   A simpler alternative to using `MatRegister()` for an application specific matrix format is to use `MatCreateShell()`, which
    280:   generates a `Mat` of `MatType` `MATSHELL`. One can then use `MatShellSetContext()` and `MatShellSetOperation()` to provide
    281:   the data structures and routines customized for their matrix format.

    283:   Example Usage:
    284: .vb
    285:    MatRegister("my_mat", MyMatCreate);
    286: .ve

    288:   Then, your solver can be chosen with the procedural interface via `MatSetType(Mat, "my_mat")` or at runtime via the option
    289:   `-mat_type my_mat`

    291: .seealso: [](ch_matrices), `Mat`, `MatType`, `MatSetType()`, `MatRegisterAll()`
    292: @*/
    293: PetscErrorCode MatRegister(const char sname[], PetscErrorCode (*function)(Mat))
    294: {
    295:   PetscFunctionBegin;
    296:   PetscCall(MatInitializePackage());
    297:   PetscCall(PetscFunctionListAdd(&MatList, sname, function));
    298:   PetscFunctionReturn(PETSC_SUCCESS);
    299: }

    301: MatRootName MatRootNameList = NULL;

    303: /*@
    304:   MatRegisterRootName - Registers a name that can be used for either a sequential or its corresponding parallel matrix type.

    306:   Input Parameters:
    307: + rname - the rootname, for example, `MATAIJ`
    308: . sname - the name of the sequential matrix type, for example, `MATSEQAIJ`
    309: - mname - the name of the parallel matrix type, for example, `MATMPIAIJ`

    311:   Level: developer

    313:   Notes:
    314:   `MatSetType()` and `-mat_type name` will automatically use the sequential or parallel version
    315:   based on the size of the MPI communicator associated with the matrix.

    317:   The matrix rootname should not be confused with the base type of the function
    318:   `PetscObjectBaseTypeCompare()`

    320:   Developer Notes:
    321:   PETSc vectors have a similar rootname that indicates PETSc should automatically select the appropriate `VecType` based on the
    322:   size of the communicator but it is implemented by simply having additional `VecCreate_RootName()` registerer routines that dispatch to the
    323:   appropriate creation routine. Why have two different ways of implementing the same functionality for different types of objects? It is
    324:   confusing.

    326: .seealso: [](ch_matrices), `Mat`, `MatType`, `PetscObjectBaseTypeCompare()`
    327: @*/
    328: PetscErrorCode MatRegisterRootName(const char rname[], const char sname[], const char mname[])
    329: {
    330:   MatRootName names;

    332:   PetscFunctionBegin;
    333:   PetscCall(PetscNew(&names));
    334:   PetscCall(PetscStrallocpy(rname, &names->rname));
    335:   PetscCall(PetscStrallocpy(sname, &names->sname));
    336:   PetscCall(PetscStrallocpy(mname, &names->mname));
    337:   if (!MatRootNameList) {
    338:     MatRootNameList = names;
    339:   } else {
    340:     MatRootName next = MatRootNameList;
    341:     while (next->next) next = next->next;
    342:     next->next = names;
    343:   }
    344:   PetscFunctionReturn(PETSC_SUCCESS);
    345: }